#ifndef _QGCM_DIAGFUNCS_SRC_
#define _QGCM_DIAGFUNCS_SRC_
#include "diagfuncs.h"
#include "qgcmio.h"
#include "main_aux.h"
#include <stdarg.h>
#include <mpValReader.h>
#include <mpParser.h>
#include <limits>
#include <signal.h>


using namespace mup;
// Globals ** :

handle_collection handler;
qstream stored_diagout;
/* ------------------------------------------------------ */
/* Global scope functions: */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
int mpError(const char * format, ...){
  // variadic throws ParserError e with printf formatting...
  // typically the model will terminate after this
  va_list(args);
  va_start(args, format);
  string msg;
  int num = stdvsprintf(msg, format, args);
  throw ParserError(msg);
  return num;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
string & local_get_errmsg_str(){
  static string stored_msg;
  return stored_msg;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void set_mp_errmsg(const char * fmt, ...){
  va_list(args);
  va_start(args, fmt);
  string & msg = local_get_errmsg_str();  
  stdsprintf(msg, fmt, args);
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
const string & get_mp_errmsg(){
  return local_get_errmsg_str();
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void set_default_value_for_type(IValue & val){
  const IValue & def = get_default_value_of_type(val.GetType());
  val = def;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
string gethandleinfo( const IValue & val ){
  if(is_value_valid_handle(val)){
    const qstream * fout = handler.get_file( val );
    if(fout){
      return "(FILE-->\"" + fout->getpath() + string("\")");
    }else{
      return "(null)";
    }
  }else{
    return "(invalid)";
  }
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_basic_parser(ParserX & p, const string & type_list){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_bool = ( type_list.find('b') != string::npos )|| include_all;
  bool include_string = ( type_list.find('s') != string::npos )|| include_all;
  bool include_complex = ( type_list.find('c') != string::npos ) || include_all;
  // define basic parser stuff
  p.AddValueReader(new DblValReader());
  if(include_bool) p.AddValueReader(new BoolValReader());
  if(include_string) p.AddValueReader(new StrValReader());
  // Constants
  p.DefineConst( _T("Pi"), static_cast<float_type>( 3.141592653589793238462643 ) );
  p.DefineConst( _T("PI"), static_cast<float_type>( 3.141592653589793238462643 ) );
  p.DefineConst( _T("e"),  static_cast<float_type>( 2.718281828459045235360287 ) );
  if( include_complex) p.DefineConst( _T("I"), cmplx_type(0, 1) ); 
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_self_assign_operators(model_class * mo, ParserX & p, const string & type_list, bool optimized){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_int = ( type_list.find('i') != string::npos ) || include_all;

  if(include_int){ 
    pair<string, diag_self_assign_oprt_func> cstack( "i:i", diag_oprt_self_assign_increment_int);
    p.DefineInfixOprt( new diagSelfAssignOprt(mo, "++", cstack, optimized) );
  }

  if(include_int) { 
    //flist.push_back( pair<string, diag_self_assign_oprt_func>( "i:i", diag_oprt_self_assign_decrement_int) );
    pair<string, diag_self_assign_oprt_func> cstack( "i:i", diag_oprt_self_assign_decrement_int);
    p.DefineInfixOprt( new diagSelfAssignOprt(mo, "--", cstack, optimized) ); 
  }
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assignBuiltIns(ParserX & p, model_class * mo, const char * name, vector< pair<string, diag_builtin_func> > & flist, bool optimized)
{
  for(int i = 0; i < flist.size(); i++){
    p.DefineFun( new diagBuiltInFunc(mo, name, flist[i], optimized) );
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_builtin_functions(model_class * mo, ParserX & p, const string & type_list, bool optimized){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_int = ( type_list.find('i') != string::npos ) || include_all;
  bool include_bool = ( type_list.find('b') != string::npos )|| include_all;
  bool include_real = ( type_list.find('f') != string::npos )|| include_all;
  bool include_complex = ( type_list.find('c') != string::npos )|| include_all;
  bool include_string = ( type_list.find('s') != string::npos )|| include_all;
  bool include_array = ( type_list.find('a') != string::npos )|| include_all;

  pair< string, diag_builtin_func> dpair;
  vector< pair<string, diag_builtin_func> > flist;
#define add_func(list,func) { dpair.first = list; dpair.second=func; flist.push_back(dpair); };
  // abs 
  if(include_int)     add_func("i:i", diag_builtin_abs_int);
  if(include_real)    add_func("f:f", diag_builtin_abs_real);
  if(include_complex) add_func("c:c", diag_builtin_abs_complex);
  assignBuiltIns(p, mo, "abs", flist, optimized);
  flist.clear();
  // sin
  if(include_real)    add_func("f:f", diag_builtin_sin_real);
  (assignBuiltIns(p, mo, "sin", flist, optimized) );
  flist.clear();
  // cos
  if(include_real)    add_func("f:f", diag_builtin_cos_real);
  (assignBuiltIns(p, mo, "cos", flist, optimized) );
  flist.clear();
  // tan
  if(include_real)    add_func("f:f", diag_builtin_tan_real);
  (assignBuiltIns(p, mo, "tan", flist, optimized) );
  flist.clear();
  // asin
  if(include_real)    add_func("f:f", diag_builtin_asin_real);
  (assignBuiltIns(p, mo, "asin", flist, optimized) );
  flist.clear();
  // acos
  if(include_real)    add_func("f:f", diag_builtin_acos_real);
  (assignBuiltIns(p, mo, "acos", flist, optimized) );
  flist.clear();
  // atan
  if(include_real)    add_func("f:f", diag_builtin_atan_real);
  (assignBuiltIns(p, mo, "atan", flist, optimized) );
  flist.clear();
  // sinh
  if(include_real)    add_func("f:f", diag_builtin_sinh_real);
  (assignBuiltIns(p, mo, "sinh", flist, optimized) );
  flist.clear();
  // cosh
  if(include_real)    add_func("f:f", diag_builtin_cosh_real);
  (assignBuiltIns(p, mo, "cosh", flist, optimized) );
  flist.clear();
  // tanh
  if(include_real)    add_func("f:f", diag_builtin_tanh_real);
  (assignBuiltIns(p, mo, "tanh", flist, optimized) );
  flist.clear();
  // natural log
  if(include_real)    add_func("f:f", diag_builtin_ln_real);
  (assignBuiltIns(p, mo, "log", flist, optimized) );
  (assignBuiltIns(p, mo, "ln", flist, optimized) );
  flist.clear();
  // log10
  if(include_real)    add_func("f:f", diag_builtin_log10_real);
  (assignBuiltIns(p, mo, "log10", flist, optimized) );
  flist.clear();
  // exp
  if(include_real)    add_func("f:f", diag_builtin_exp_real);
  if(include_complex) add_func("c:c", diag_builtin_exp_complex);
  (assignBuiltIns(p, mo, "exp", flist, optimized) );
  flist.clear();
  // sqrt
  if(include_complex) add_func("f:f", diag_builtin_sqrt_complex);
  if(include_real) add_func("f:f", diag_builtin_sqrt_real);
  (assignBuiltIns(p, mo, "sqrt", flist, optimized) );
  flist.clear();
  // re
  if(include_complex) add_func("f:c", diag_builtin_real_complex);
  (assignBuiltIns(p, mo, "Re", flist, optimized) );
  flist.clear();
  // Im
  if(include_complex) add_func("f:c", diag_builtin_imag_complex);
  (assignBuiltIns(p, mo, "Im", flist, optimized) );
  flist.clear();
  // Conj
  if(include_complex) add_func("c:c", diag_builtin_conj_complex);
  (assignBuiltIns(p, mo, "conj", flist, optimized) );
  flist.clear();
  // min
  if(include_int)  add_func("i:ii*", diag_builtin_min_int);
  if(include_real) add_func("f:ff*", diag_builtin_min_real);
  (assignBuiltIns(p, mo, "min", flist, optimized) );
  flist.clear();
  // max
  if(include_int)  add_func("i:ii*", diag_builtin_max_int);
  if(include_real) add_func("f:ff*", diag_builtin_max_real);
  (assignBuiltIns(p, mo, "max", flist, optimized) );
  flist.clear();
  // sum
  if(include_int)  add_func("i:ii*", diag_builtin_sum_int);
  if(include_real) add_func("f:ff*", diag_builtin_sum_real);
  (assignBuiltIns(p, mo, "sum", flist, optimized) );
  flist.clear();
  // length
  if(include_string)  add_func("i:s", diag_builtin_length_string);
  if(include_array)   add_func("i:a", diag_builtin_length_array);
  (assignBuiltIns(p, mo, "length", flist, optimized) );
  flist.clear();

  if(include_array)   add_func("i:a", diag_builtin_length_array);
  (assignBuiltIns(p, mo, "size", flist, optimized) );

#undef add_func
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assignDiagOprtBin(ParserX & p, model_class * mo, const char * name, 
    EOprtPrecedence priority, EOprtAsct assoc,
    vector< pair<string, diag_bin_oprt_func> > & flist, bool optimized)
{
  for(int i = 0; i < flist.size(); i++){
    p.DefineOprt( new diagOprtBin(mo, name, priority, assoc, flist[i], optimized) );
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_binary_operators(model_class * mo, ParserX & p, const string & type_list, bool optimized){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_int = ( type_list.find('i') != string::npos ) || include_all;
  bool include_all_int = ( type_list.find('i') != string::npos ) || include_all;
  bool include_bool = ( type_list.find('b') != string::npos )|| include_all;
  bool include_real = ( type_list.find('f') != string::npos )|| include_all;
  bool include_complex = ( type_list.find('c') != string::npos )|| include_all;
  bool include_string = ( type_list.find('s') != string::npos )|| include_all;
  bool include_array = ( type_list.find('a') != string::npos )|| include_all;

  vector< pair< string, diag_bin_oprt_func > > flist;
#define add_func(cstack,func) { flist.push_back( pair<string,diag_bin_oprt_func>(string(cstack),func) ); assert(strlen(cstack)==4); }
  // operator+
  if(include_int)       add_func("i:ii", diag_bin_oprt_add_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_add_real);
  if(include_complex)   add_func("c:cc", diag_bin_oprt_add_complex);
  if(include_string)    add_func("s:ss", diag_bin_oprt_add_string);
  (assignDiagOprtBin(p, mo,"+",prADD_SUB,oaLEFT,flist,optimized) );
  flist.clear();
  // operator-
  if(include_int)       add_func("i:ii", diag_bin_oprt_sub_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_sub_real);
  if(include_complex)   add_func("c:cc", diag_bin_oprt_sub_complex);
  (assignDiagOprtBin(p, mo,"-",prADD_SUB,oaLEFT,flist,optimized) );
  flist.clear();
  // operator*
  if(include_all_int)   add_func("i:ii", diag_bin_oprt_mul_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_mul_real);
  if(include_complex)   add_func("c:cc", diag_bin_oprt_mul_complex);
  (assignDiagOprtBin(p, mo,"*",prMUL_DIV,oaLEFT,flist,optimized) );
  flist.clear();
  // operator/
  if(include_all_int)   add_func("i:ii", diag_bin_oprt_div_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_div_real);
  if(include_complex)   add_func("c:cc", diag_bin_oprt_div_complex);
  (assignDiagOprtBin(p, mo,"/",prMUL_DIV,oaLEFT,flist,optimized) );
  flist.clear();
  // operator%
  if(include_int)      add_func("i:ii", diag_bin_oprt_mod_int);
  if(include_real && ! include_int ) add_func("f:ff", diag_bin_oprt_mod_real);
  (assignDiagOprtBin(p, mo,"%",prMUL_DIV,oaLEFT,flist,optimized) );
  flist.clear();
  // operator^
  //if(include_all_int)   add_func("i:ii", diag_bin_oprt_pow_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_pow_real);
  if(include_complex)   add_func("c:cc", diag_bin_oprt_pow_complex);
  (assignDiagOprtBin(p, mo,"^",prPOW,oaRIGHT,flist,optimized) );
  flist.clear();


  // operator #
  if(include_int && include_array) add_func("f:ai", diag_bin_oprt_sharp);
  // has higher priority than all binary operators
  (assignDiagOprtBin(p, mo,"#",static_cast<EOprtPrecedence>(11),oaLEFT,flist,optimized) );
  flist.clear();

  // ops that return type bool but take other vals:
  if(include_bool){
    // operator==
    add_func("b:bb",diag_bin_oprt_eq_bool);
    if(include_int)     add_func("b:ii",diag_bin_oprt_eq_int);
    if(include_real)    add_func("b:ff",diag_bin_oprt_eq_real);
    if(include_complex) add_func("b:cc",diag_bin_oprt_eq_complex);
    if(include_string)  add_func("b:ss",diag_bin_oprt_eq_string);
    if(include_array)   add_func("b:aa",diag_bin_oprt_eq_array);
    (assignDiagOprtBin(p, mo,"==",prRELATIONAL1,oaLEFT,flist,optimized) );
    flist.clear();
    // operator!= 
    add_func("b:bb",diag_bin_oprt_neq_bool);
    if(include_int)     add_func("b:ii",diag_bin_oprt_neq_int);
    if(include_real)    add_func("b:ff",diag_bin_oprt_neq_real);
    if(include_complex) add_func("b:cc",diag_bin_oprt_neq_complex);
    if(include_string)  add_func("b:ss",diag_bin_oprt_neq_string);
    if(include_array)   add_func("b:aa",diag_bin_oprt_neq_array);
    (assignDiagOprtBin(p, mo,"!=",prRELATIONAL1,oaLEFT,flist,optimized) );
    flist.clear();
    /// operator<=
    if(include_int)     add_func("b:ii",diag_bin_oprt_leq_int);
    if(include_real)    add_func("b:ff",diag_bin_oprt_leq_real);
    (assignDiagOprtBin(p, mo,"<=",prRELATIONAL2,oaLEFT, flist, optimized) );
    (assignDiagOprtBin(p, mo,"=<",prRELATIONAL2,oaLEFT, flist, optimized) );
    flist.clear();
    /// operator=>
    if(include_int)     add_func("b:ii",diag_bin_oprt_geq_int);
    if(include_real)    add_func("b:ff",diag_bin_oprt_geq_real);
    (assignDiagOprtBin(p, mo,"=>",prRELATIONAL2,oaLEFT, flist, optimized) );
    (assignDiagOprtBin(p, mo,">=",prRELATIONAL2,oaLEFT, flist, optimized) );
    flist.clear();
    /// operator<
    if(include_int)     add_func("b:ii",diag_bin_oprt_lt_int);
    if(include_real)    add_func("b:ff",diag_bin_oprt_lt_real);
    (assignDiagOprtBin(p, mo,"<",prRELATIONAL2,oaLEFT, flist, optimized) );
    flist.clear();
    /// operator>
    if(include_int)     add_func("b:ii",diag_bin_oprt_gt_int);
    if(include_real)    add_func("b:ff",diag_bin_oprt_gt_real);
    (assignDiagOprtBin(p, mo,">",prRELATIONAL2,oaLEFT, flist, optimized) );
    flist.clear();
    // boolean operator &&
    add_func("b:bb",diag_bin_oprt_and_bool);
    (assignDiagOprtBin(p, mo,"&&",prLOGIC_AND,oaLEFT,flist,optimized) );
    flist.clear();
    // boolean operator || 
    add_func("b:bb",diag_bin_oprt_or_bool);
    (assignDiagOprtBin(p, mo,"||",prLOGIC_OR,oaLEFT,flist,optimized) );
    flist.clear();
  }

#undef add_func
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assignprofileOprtBin(ParserX & p, model_class * mo, const char * name, 
    EOprtPrecedence priority, EOprtAsct assoc,
    vector< pair< string, diag_bin_oprt_func> > & flist, const int pindex){
  for(int i = 0; i < flist.size(); i++){
    p.DefineOprt( new profileOprtBin(mo, name, priority, assoc, flist[i], pindex) );
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_profile_binary_operators(model_class * mo, ParserX & p, const string & type_list, const int pindex){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_int = ( type_list.find('i') != string::npos ) || include_all;
  bool include_all_int = ( type_list.find('i') != string::npos );
  // bool include_bool = ( type_list.find('b') != string::npos )|| include_all; (no boolean support right now)
  bool include_real = ( type_list.find('f') != string::npos )|| include_all;
  bool include_string = ( type_list.find('s') != string::npos )|| include_all;
  // bool include_array = ( type_list.find('a') != string::npos )|| include_all; // no direct array manipulation

  vector< pair< string, diag_bin_oprt_func > > flist;
#define add_func(cstack,func) { flist.push_back( pair<string,diag_bin_oprt_func>(string(cstack),func) ); assert(strlen(cstack)==4); }
  // operator+
  //if(include_int)       add_func("i:ii", diag_bin_oprt_add_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_add_real);
  if(include_real)      add_func("f:fi", diag_bin_oprt_add_real);
  if(include_real)      add_func("f:if", diag_bin_oprt_add_real);
  if(include_real)      add_func("f:af", diag_bin_oprt_add_real);
  if(include_real)      add_func("f:fa", diag_bin_oprt_add_real);
  if(include_string)    add_func("s:ss", diag_bin_oprt_add_string);
  if(include_string)    add_func("s:sa", diag_bin_oprt_add_string);
  if(include_string)    add_func("s:as", diag_bin_oprt_add_string);
  (assignprofileOprtBin(p, mo,"+",prADD_SUB,oaLEFT,flist,pindex) );
  flist.clear();
  // operator-
  //if(include_int)       add_func("i:ii", diag_bin_oprt_sub_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_sub_real);
  if(include_real)      add_func("f:fa", diag_bin_oprt_sub_real);
  if(include_real)      add_func("f:af", diag_bin_oprt_sub_real);
  (assignprofileOprtBin(p, mo,"-",prADD_SUB,oaLEFT,flist,pindex) );
  flist.clear();
  // operator*
  //if(include_all_int)   add_func("i:ii", diag_bin_oprt_mul_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_mul_real);
  if(include_real)      add_func("f:fi", diag_bin_oprt_mul_real);
  if(include_real)      add_func("f:if", diag_bin_oprt_mul_real);
  if(include_real)      add_func("f:fa", diag_bin_oprt_mul_real);
  if(include_real)      add_func("f:af", diag_bin_oprt_mul_real);
  (assignprofileOprtBin(p, mo,"*",prMUL_DIV,oaLEFT,flist,pindex) );
  flist.clear();
  // operator/
  //if(include_all_int)   add_func("i:ii", diag_bin_oprt_div_int);
  if(include_real)      add_func("f:ff", diag_bin_oprt_div_real);
  if(include_real)      add_func("f:fi", diag_bin_oprt_div_real);
  if(include_real)      add_func("f:if", diag_bin_oprt_div_real);
  if(include_real)      add_func("f:af", diag_bin_oprt_div_real);
  if(include_real)      add_func("f:fa", diag_bin_oprt_div_real);
  (assignprofileOprtBin(p, mo,"/",prMUL_DIV,oaLEFT,flist,pindex) );
  flist.clear();
  // operator%
  if(include_int)      add_func("i:ii", diag_bin_oprt_mod_int);
  if(include_int)      add_func("i:ia", diag_bin_oprt_mod_int);
  if(include_int)      add_func("i:ai", diag_bin_oprt_mod_int);
  (assignprofileOprtBin(p, mo,"%",prMUL_DIV,oaLEFT,flist,pindex) );
  flist.clear();

  // operator^

  if(include_real)      add_func("f:ff", diag_bin_oprt_pow_real);
  if(include_real)      add_func("f:fi", diag_bin_oprt_pow_real);
  if(include_real)      add_func("f:if", diag_bin_oprt_pow_real);
  if(include_real)      add_func("f:fa", diag_bin_oprt_pow_real);
  if(include_real)      add_func("f:af", diag_bin_oprt_pow_real);
  (assignprofileOprtBin(p, mo,"^",prPOW,oaRIGHT,flist,pindex) );
  flist.clear();


  // operator #
  add_func("i:ai", profile_bin_oprt_sharp);
  // has higher priority than all binary operators
  (assignprofileOprtBin(p, mo,"#",static_cast<EOprtPrecedence>(11),oaLEFT,flist, -1) );
  flist.clear();


#undef add_func
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assigndiagAssignOprt(ParserX & p, model_class * mo, const char * name, 
    EOprtPrecedence priority, EOprtAsct assoc,
    vector< pair< string, diag_assign_oprt_func> > & flist, bool optimized){
  for(int i = 0; i < flist.size(); i++){
    p.DefineOprt( new diagAssignOprt(mo, name, priority, assoc, flist[i], optimized) );
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_assign_operators(model_class * mo, ParserX & p, const string & type_list, bool optimized){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_int = ( type_list.find('i') != string::npos ) || include_all;
  bool include_bool = ( type_list.find('b') != string::npos )|| include_all;
  bool include_real = ( type_list.find('f') != string::npos )|| include_all;
  bool include_complex = ( type_list.find('c') != string::npos )|| include_all;
  bool include_string = ( type_list.find('s') != string::npos )|| include_all;
  bool include_array = ( type_list.find('a') != string::npos )|| include_all;
  //bool include_handle = ( type_list.find('h') != string::npos ) || include_all;

  vector< pair< string, diag_assign_oprt_func > > flist;
#define add_func(cstack,func) { flist.push_back( pair<string,diag_assign_oprt_func>(string(cstack),func) ); assert(strlen(cstack)==4); }
  // operator=
  if(include_bool)    add_func("b:bb",diag_oprt_assign_bool_2_bool);
  if(include_int)     add_func("i:ii",diag_oprt_assign_int_2_int);
  if(include_real)    add_func("f:fi",diag_oprt_assign_real_2_real);
  if(include_real)    add_func("f:ff",diag_oprt_assign_real_2_real);
  if(include_complex) add_func("c:cc",diag_oprt_assign_complex_2_complex);
  if(include_string)  add_func("s:ss",diag_oprt_assign_string_2_string);
  if(include_array)   add_func("a:aa",diag_oprt_assign_array_2_array);
  if(include_all)     add_func("v:vv",diag_oprt_assign_variable_2_variable); // typeless assignment
  (assigndiagAssignOprt(p, mo, "=", prASSIGN, oaLEFT, flist, optimized) );
  flist.clear();
  // operator +=
  if(include_int)     add_func("i:ii",diag_oprt_assign_add_int);
  if(include_real)    add_func("f:ff",diag_oprt_assign_add_real);
  if(include_complex) add_func("c:cc",diag_oprt_assign_add_complex);
  if(include_string)  add_func("s:ss",diag_oprt_assign_add_string);
  (assigndiagAssignOprt(p, mo, "+=", prASSIGN, oaLEFT, flist, optimized) );
  flist.clear(); 
  // operator -=
  if(include_int)     add_func("i:ii",diag_oprt_assign_sub_int);
  if(include_real)    add_func("f:ff",diag_oprt_assign_sub_real);
  if(include_complex) add_func("c:cc",diag_oprt_assign_sub_complex);
  (assigndiagAssignOprt(p, mo, "-=", prASSIGN, oaLEFT, flist, optimized) );
  flist.clear();
  // operator *=
  if(include_int)     add_func("i:ii",diag_oprt_assign_mul_int);
  if(include_real)    add_func("f:ff",diag_oprt_assign_mul_real);
  if(include_complex) add_func("c:cc",diag_oprt_assign_mul_complex);
  (assigndiagAssignOprt(p, mo, "*=", prASSIGN, oaLEFT, flist, optimized) );
  flist.clear();
  // operator /=
  if(include_int)     add_func("i:ii",diag_oprt_assign_div_int);
  if(include_real)    add_func("f:ff",diag_oprt_assign_div_real);
  if(include_complex) add_func("c:cc",diag_oprt_assign_div_complex);
  (assigndiagAssignOprt(p, mo, "/=", prASSIGN, oaLEFT, flist, optimized) );
  flist.clear();
#undef add_func
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assigndiagInfixOprt(ParserX & p, model_class * mo, const char * name, 
    vector< pair< string, diag_infix_oprt_func> > & flist, bool optimized){
  for(int i = 0; i < flist.size(); i++){
    p.DefineInfixOprt( new diagInfixOprt(mo, name, flist[i], optimized) );
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void include_infix_operators(model_class * mo, ParserX & p, const string & type_list, bool optimized){
  bool include_all = ( type_list.find('v') != string::npos );
  bool include_int = ( type_list.find('i') != string::npos ) || include_all;
  bool include_bool = ( type_list.find('b') != string::npos )|| include_all;
  bool include_real = ( type_list.find('f') != string::npos )|| include_all;
  bool include_complex = ( type_list.find('c') != string::npos )|| include_all;
  bool include_string = ( type_list.find('s') != string::npos )|| include_all;
  bool include_array = ( type_list.find('a') != string::npos )|| include_all;

  vector< pair< string, diag_infix_oprt_func> > flist;
  if(include_bool) flist.push_back( pair<string, diag_infix_oprt_func>( "b:b", diag_uni_oprt_not_bool ) );
  (assigndiagInfixOprt(p, mo, "!", flist, optimized ) );

  flist.clear();
  if( include_int ) flist.push_back( pair<string, diag_infix_oprt_func>( "i:i", diag_uni_oprt_neg_int ) );
  if( include_real ) flist.push_back( pair<string, diag_infix_oprt_func>( "f:f", diag_uni_oprt_neg_float ) );
  if( include_complex ) flist.push_back( pair<string, diag_infix_oprt_func>( "c:c", diag_uni_oprt_neg_complex ) );
  (assigndiagInfixOprt(p, mo, "-", flist, optimized) );

};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void parser_define_basics(ParserX & p, model_class * mo, const char * types, bool optimized, bool diag, int pindex){
  string stypes(types);
  rassert( stypes.size() > 0 );
  bool include_all = ( stypes.find('v') != string::npos );
  bool include_array = ( stypes.find('a') != string::npos )|| include_all;
  include_basic_parser(p, stypes );
  include_assign_operators(mo, p, stypes, optimized );
  include_self_assign_operators(mo, p, stypes, optimized );
  include_infix_operators(mo, p, stypes, optimized );
  if( diag ){
    include_binary_operators( mo, p , stypes, optimized );
  }else{
    include_profile_binary_operators(mo, p, stypes, pindex);
  }
  include_builtin_functions(mo, p, stypes, optimized );
  // add the array indexing operator function myarray[index]
  if(include_array)  p.DefineArrayOprt( diagOprtArrayIndex ); 
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
const IValue & get_default_value_of_type(char_type ct){
  static char_type handle_type = 'h';
  static Value int_val(static_cast<int_type>(1));
  static Value flt_val(static_cast<float_type>(0.1));
  static Value cp_val( cmplx_type(1,1) );
  static Value str_val( string_type("") );
  static Value handle_val( handle_type );
  static Value bool_val( static_cast<bool_type>(false) );
  static Value var_val; // set to 'v' 
  static Value array_val( vector<Value>(1, flt_val) );
  assert( var_val.GetType() == 'v' );
  switch( ct ){
    case 'i' : return * int_val.AsIValue();
    case 'f' : return * flt_val.AsIValue();
    case 'c' : return * cp_val.AsIValue();
    case 's' : return * str_val.AsIValue();
    case 'h' : return * handle_val.AsIValue();
    case 'b' : return * bool_val.AsIValue();
    case 'a' : return * array_val.AsIValue();
    case 'v' : return * var_val.AsIValue();
    default: return * var_val.AsIValue();
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool is_value_valid_handle(const IValue & val){
  if( val.GetType() != 'i' &&  val.GetType() != 'h' ) return false;
  int ival = static_cast<int>( val.GetInteger() );
  if(ival == 0) return true; // null is a valid handle
  // check the handler stack
  const qstream * fout = handler.get_file( val );
  if(fout==NULL) return false;
  // return true only if the handle is attached to an open file
  return ! (fout->is_closed());
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool is_type_valid(const char_type ct){
  return (ct == 'b' || ct=='a'||ct == 'c'|| ct == 'i' || ct == 'f' || ct == 's' || ct == 'h' || ct == 'v' );
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void create_function_call_stack_string(const vector<char_type> & stacklist, string & dest){
	dest.push_back('(');
	for(int i = 0; i < static_cast<int>(stacklist.size()); i++){
		dest+= get_type_id( stacklist[i] );
		if( i < static_cast<int>(stacklist.size()) - 1 ) dest+=", ";
	}
	dest.push_back(')');
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
char_type get_var_type_from_string(const string & str){
  if( str == "array" ) return 'a';
  if( str == "complex" ) return 'c';
  if( str == "real" ) return 'f';
  if( str == "string" ) return 's';
  if( str == "handle" ) return 'h';
  if( str == "variable" ) return 'v';
  if( str == "integer" ) return 'i';
  if( str == "bool" ) return 'b';
  mpError("%s does not name a type.", str.c_str());
  return 0;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool is_string_var_type(const string & str){
  if( str == "array" ) return true;
  if( str == "complex" ) return true;
  if( str == "integer" ) return true;
  if( str == "real" ) return true;
  if( str == "string" ) return true;
  if( str == "handle" ) return true;
  if( str == "bool" ) return true;
  if( str == "variable" ) return true;
  return false;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool does_type_conform_to_type(char_type t1, char_type t2)
{
  // does t1 conform to t2?
  if( t2 == 'v' ) return true;
  switch(t2){
    case 'h':
      return t1 == 'h';
    case 'i':
      return t1 == 'i' || t1 == 'b';
    case 'f':
      return t1 == 'i' || t1 == 'b' || t1 == 'f';
    case 'c':
      return t1 == 'f' || t1 == 'i' || t1 == 'b' || t1 == 'c';
    case 'b':
      return t1 == 'b';
    case 's':
      return t1 == 's';
    case 'a':
      return t1 == 'a';
    default :
      rassert("Unrecognized type.");
      return false;
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool does_value_conform_to_type(const IValue & val, char_type t){
  if(t == 'v') return true;
  switch(t){
    case 'h':
      return is_value_valid_handle(val);
    case 'i':
      return val.GetType() == 'i' && ((val.GetFloat() != val.GetFloat()) || val.IsInteger() || val.GetType() == 'b'); 
      // allow Nan cases
    case 'f':
      return (val.GetFloat() != val.GetFloat()) || val.IsInteger() || val.GetType() == 'f' || val.GetType() == 'b';
    case 'c':
      return val.GetType() == 'c' || val.GetType() == 'i' || val.GetType() == 'f' || val.GetType() == 'b';
    case 'b':
      return val.GetType() == 'b';
    case 's':
      return val.GetType() == 's';
    case 'a':
      return val.GetType() == 'a';
    default:
      return false; 
  }
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
char cast_to_char(const IValue & val){
  char_type type = val.GetType();
  if(type == 'i'){
    return static_cast<char>(val.GetInteger());
  }else if(type == 'f'){
    return static_cast<char>(val.GetFloat());
  }else if(type == 's'){
    return static_cast<char>(convertstring2integer(val.GetString()));
  }else if(type == 'b'){
    return val.GetBool() ? 1 : 0;
  }else{
    // conversion not supported
    mpError("Cannot convert type %s to char.", get_type_id(type));
    return 0;
  }
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
int cast_to_int(const IValue & val){
  char_type type = val.GetType();
  if(type == 'i'){
    return static_cast<int>(val.GetInteger());
  }else if(type == 'f'){
    return static_cast<int>(val.GetFloat());
  }else if(type == 's'){
    return static_cast<int>(convertstring2integer(val.GetString()));
  }else if(type == 'b'){
    return val.GetBool() ? 1 : 0;
  }else{
    // conversion not supported
    mpError("Cannot convert type %s to int.", get_type_id(type));
    return 0;
  }
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
float cast_to_float(const IValue & val){
  char_type type = val.GetType();
  if(type == 'i'){
    return static_cast<float>(val.GetInteger());
  }else if(type == 'c'){
    return static_cast<float>(val.GetComplex().real());
  }else if(type == 'f'){
    return static_cast<float>(val.GetFloat());
  }else if(type == 's'){
    return static_cast<float>(convertstring2float(val.GetString()));
  }else if(type == 'b'){
    return val.GetBool() ? 1.0 : 0.0;
  }else{
    // conversion not supported
    mpError("Cannot convert type %s to real.", get_type_id(type));
    return 0;
  }
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
string cast_to_string(const IValue & val){
  // all conversions are supported
  char_type type = val.GetType();
  //char * bufr;
  if(type == 'v') return string("<empty>");
  if(type == 'a') return string("array(size=") + convertint2string(val.GetArray().size()) + string(")");
  if(type == 'i') return convertint2string(static_cast<int>(val.GetInteger()));
  if(type == 'f'){
    return convertfloat2string( static_cast<float>( val.GetFloat() ) );
  }
  if(type == 'c'){
    const cmplx_type & cp = val.GetComplex();
    return convertcomplex2string( qgcomplex(cp.real(), cp.imag() ) );
  }
  if(type == 'h'){
    return converthandle2string( val.GetInteger() );
  }
  if(type == 's') return val.GetString();
  if(type == 'b') {
    if(val.GetBool()) return string("true");
    return string("false");
  }
  return string("");
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool cast_to_bool(const IValue & val){
  char_type type = val.GetType();
  switch( type ){
    case 'i':
      return val.GetInteger() != 0;
    case 'f':
      return val.GetFloat() != 0;
    case 'b':
      return val.GetBool();
    default:
      mpError("Cannot convert type %s to bool.", get_type_id(type));
  }
  return false;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assert_stack_is_valid(const string & stack){
  // stack strings must have the form (return type):(argt1, argt2, ...)
  // i.e. i:ffc or i:sv* 
  // * means variadic (no type matching can be done);
  rassert( stack.size() >= 2 );
  rassert( stack[1] == ':' );
  rassert( is_type_valid( stack[0] ) ); // the return type
  for(size_t iter = 2; iter < stack.size(); iter++){
    if( stack[iter] == '*' ){
      rassert( iter == stack.size() - 1 );
    }else{
      rassert( is_type_valid( stack[iter] ) );
    }
  }
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
int get_stack_argc(const string & stack){
  assert_stack_is_valid(stack);
  if( stack[stack.size()-1]  == '*' ) return -1;
  return static_cast<int>(stack.size()) - 2;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void assert_stack_can_be_optimized(const string & cstack){
  char_type rtype = 0;
  for(size_t iter = 0; iter < cstack.size(); iter++){
    rtype = cstack[iter];
    if(iter != 1){ 
      if( ! ( rtype=='c' || rtype=='i' || rtype=='f' ) ){ 
        fatal_stop("Call stack cannot be optimized.");
      }
    }
  }
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool does_stack_match(const string & fcstack, const ptr_val_type * argv, const int argc){
#ifdef DEBUG
  assert_stack_is_valid(fcstack);
#endif
  const int clen = fcstack.size() - 2; // ignore the return and colon, ex: i:fffa -> fffa
  const char * cstack = &((fcstack.c_str())[2]);
  for(int iter = 0; iter < clen ; iter++){
    if( cstack[iter] == '*' ) break;
    if(iter >= argc) return false;
    if( ! does_value_conform_to_type(* (argv[iter].Get()) , cstack[iter]) ) return false;
  }
  return true;
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
bool does_return_match(const string & fcstack, const ptr_val_type & ret){
#ifdef DEBUG
  assert_stack_is_valid(fcstack);
#endif
  return does_value_conform_to_type( *(ret.Get()), fcstack[0] );
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void gen_call_desc(string & desc, const string & type, const string & name, const string & stack){
#ifdef DEBUG
  assert_stack_is_valid(stack); 
#endif
  desc.clear();
  string ret_type( get_type_id( stack[0] ) );
  stdsprintf(desc, "%s %s %s(", type.c_str(), ret_type.c_str(), name.c_str()); 
  for(size_t iter = 2; iter < stack.size(); iter++ ){
    if(iter > 2) desc+=string(", ");
    if( stack[iter] == '*' ){ desc += string("..."); break; }
    desc+= string(get_type_id(stack[iter]));
  }
  stdsprintf(desc, "%s)", desc.c_str());
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void throw_return_error(const ICallback * cb, const ptr_val_type & ret,  const ptr_val_type * argv, const int argc){
  const string & name = cb->GetIdent();
  string bad_stack;
  for(int i = 0; i < argc; i++){
    if( i > 0 ) bad_stack+=string(", ");
    bad_stack+=string(get_type_id( argv[i].Get()->GetType() ) );
  }
  mpError("\nNo matching Function: %s(%s) returns %s", name.c_str(), bad_stack.c_str(), get_type_id(ret.Get()->GetType()));
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void throw_prototype_error(const string & cbName, const string & fStack, bool isOp)
{
    rassert(false);
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void throw_prototype_error_special(const ICallback * cb, const IValue *  arg1, const IValue * arg2)
{
  const string & name = cb->GetIdent(); 
  string bad_stack;
  for(int i = 0; i < 2; i++){
    if( i > 0 ) bad_stack+=string(", ");
    bad_stack+=string(get_type_id( 
          ( i == 0 ? arg1 : arg2)->GetType() ) );
  }
  const string & possibleCalls = cb->GetDesc();
  mpError("\nNo matching operator: %s(%s)\nCandidates are: %s", 
      name.c_str(), bad_stack.c_str(), possibleCalls.c_str());
};

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
void throw_prototype_error(const ICallback * cb, const ptr_val_type * argv, const int argc, bool op=true)
{
  const string & cbName = cb->GetIdent();
  string bad_stack;
  for(int i = 0; i < argc; i++){
    if( i > 0 ) bad_stack+=string(", ");
    bad_stack+=string(get_type_id( argv[i].Get()->GetType() ) );
  }
  string type = op ? "operator" : "Function";
  const string & possibleCalls = cb->GetDesc();
  mpError("\nNo matching %s: %s(%s)\nCandidates are: %s", 
      type.c_str(), cbName.c_str(), bad_stack.c_str(), possibleCalls.c_str()
      );
};
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* : end Global scope functions */
/* ------------------------------------------------------ */

/*********************************************************************
  class handle_collection : class interface for storing of file handles 
 **********************************************************************/
handle_collection::handle_collection()
{
  int_type max_val = numeric_limits<int_type>::max();
  int_type rand_val = rand() % 1000;
  this->stdout_id = (max_val / 2) + rand_val; // starting value
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
int_type handle_collection::add_stream( qstream * new_stream ){
  // function basically decides what id to give the file
  int_type id = this->stdout_id+1;
  while( (this->file_handles.find(id) != this->file_handles.end()) ){
    id++; // increment until an empty spot is found
  }
  this->file_handles.insert( pair<int_type, qstream*>(id, new_stream) );
  return id;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Value handle_collection::file_open(const char * name, const char * fmt) {
  qstream * newstream = new qstream();
  if(newstream->open(name, fmt)){
    return Value(this->add_stream( newstream ));
  }else{
    delete newstream;
    return get_default_value_of_type('h');
  }
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const Value handle_collection::getstdout() const{
  return Value( this->stdout_id );
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
NcFile * handle_collection::get_nc_file(const IValue & handle){
  int ival = static_cast<int>(handle.GetInteger());
  return this->ncfile_handles.getncfile(ival);
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bool handle_collection::ncfile_close(const IValue & val){
  int handle = static_cast<int>(val.GetInteger());
  return this->ncfile_handles.deletencfile( handle );
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Value handle_collection::ncfile_open(const char * name, NcFile::FileMode mode){
  int_type rval = static_cast<int_type>(this->ncfile_handles.createncfile(name, mode));
  return Value(rval); // returns 0 if it doesn't work
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const qstream * handle_collection::get_file(const IValue & handle) const{
  int_type ival = handle.GetInteger();
  if( ival == stdout_id ) return &(getModel().getSTDOUT());
  map<int_type, qstream*>::const_iterator it = this->file_handles.find(ival);
  if( it == this->file_handles.end() ) return NULL;
  return it->second; 
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bool handle_collection::file_close(const IValue & handle){
  int_type ival = handle.GetInteger();
  if( ival == stdout_id ) return false;
  map<int_type, qstream*>::iterator it = this->file_handles.find(ival);
  if( it == this->file_handles.end() ) return NULL;
  qstream * fout = it->second;
  this->file_handles.erase( it );
  delete fout;
  return true;
};
/*********************************************************************
  class selfAssignOprt : class for operators of the type ++a, --a
 **********************************************************************/
diagSelfAssignOprt::diagSelfAssignOprt(model_class * model, const char * oprt_id, 
    pair< string, diag_self_assign_oprt_func> & cStack, bool optimize)
                               : IOprtInfix(oprt_id)
{
  rassert( model != NULL );
  this->isAssignmentOP = true;
  this->mo = model;
  this->func_sig = cStack.first;
  // check stack
  const int argc = get_stack_argc(func_sig);
  rassert( argc == 1 ); 
  if( optimize ) {  rassert( func_sig[2] == 'c' ); }
  this->callback = cStack.second;
  rassert( this->callback != NULL );
  this->GenDesc();
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagSelfAssignOprt::GetProtoType() const 
{
  return this->func_sig.c_str();
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagSelfAssignOprt::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  if( mo->iscompiled() ){
    this->callback(*dynamic_cast<Variable*>(argv[0].Get()) );
    return;
  }else{
    // strict type check, etc etc ...
    rassert( argc == 1 );
    const ptr_val_type & arg = argv[0];
    // first we determine which function to call based on the arg stack type
    if( does_stack_match( func_sig, argv, argc ) ){
      Variable * pVar = dynamic_cast<Variable*>( arg.Get() );
      if( ! pVar ){
        mpError("Cannot apply operator '%s' to non variable.", this->GetIdent().c_str());
      }
      callback  ( *pVar );
      *ret = *pVar;
      return;
    }
    throw_prototype_error(this, argv, argc );
  }
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagSelfAssignOprt::GenDesc(){
  string & str = this->descr;
  gen_call_desc(str, "operator", this->GetIdent(), this->func_sig);
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagSelfAssignOprt::GetDesc() const{
  return this->descr.c_str();
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken* diagSelfAssignOprt::Clone() const{
  return new diagSelfAssignOprt(*this);
}

/*********************************************************************
  class diagAssignOprt : class for operators of the type a = expr 
 **********************************************************************/
diagAssignOprt::diagAssignOprt(model_class * model, const char * token_str, EOprtPrecedence priority, EOprtAsct assoc,
    pair<string, diag_assign_oprt_func> & flist, bool optimize)
                               : IOprtBin( _T(token_str), (int)priority, assoc )
{
  rassert( model != NULL );
  this->isAssignmentOP = true;
  this->mo = model;
  this->func_sig = flist.first;
  this->callback = flist.second;
  this->optimized = optimize;
  // check stack
  int num_args = get_stack_argc( func_sig );
  rassert( num_args == 2 ); // it's a binary operator !
  if( optimized ) assert_stack_can_be_optimized( func_sig ); 
  rassert( callback != NULL );
  this->GenDesc();
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagAssignOprt::GenDesc(){
  string & str = this->descr;
  gen_call_desc(str, "operator", this->GetIdent(), this->func_sig);
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagAssignOprt::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  if( mo->iscompiled() ){
    this->callback( *dynamic_cast<Variable*>(argv[0].Get()), argv[1]  );
    return;
  }else{
    assert( argc == 2 );
    const string & cstack = func_sig;
    const ptr_val_type & arg1 = argv[0];
    const ptr_val_type & arg2 = argv[1];
    if( does_stack_match( cstack, argv, argc ) ){
      Variable *pVar = dynamic_cast<Variable*>(arg1.Get());
      if( ! pVar ){
        mpError("Cannot apply operator '%s' to non variable.", this->GetIdent().c_str());
      }
      ( callback  ) ( *pVar, arg2 );
      *ret = *pVar; // assignments also return a value which is the assigned value
      return;
    }
    throw_prototype_error(this, argv, argc );
  }
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagAssignOprt::GetProtoType() const{
  return _T(this->func_sig).c_str();
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagAssignOprt::GetDesc() const{
  return _T(this->descr.c_str());
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken* diagAssignOprt::Clone() const{
  return new diagAssignOprt(*this);
}

/*********************************************************************
  class diagInfixOprt : class interface for prefix unitary operators (i.e., !(not), -(minus) )
 **********************************************************************/
diagInfixOprt::diagInfixOprt(model_class * model, const char * oprt_id,
    pair< string, diag_infix_oprt_func> & flist, bool optimize)
                              : IOprtInfix(oprt_id)
{
  rassert( model != NULL );
  this->mo = model;
  this->func_sig = flist.first;
  this->callback = flist.second;
  // check stack
  const int argc = get_stack_argc(func_sig);
  rassert( argc == 1 );
  rassert( callback != NULL );
  if( optimize ){
    assert_stack_can_be_optimized( func_sig );    
  } 
  this->GenDesc();
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagInfixOprt::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  if( mo->iscompiled() ){
    this->callback(ret, argv[0] );
    return;
  }else{
    // strict type check, etc etc ...
    assert( argc == 1 );
    const ptr_val_type & arg = argv[0];
    // first we determine which function to call based on the arg stack type
    if( does_stack_match( func_sig, argv, argc ) ){
      ( callback ) ( ret, arg );
      return;
    } 
    // no valid operator found
    mpError("Cannot apply unitary operator '%s' to argument of type '%s'.",
        this->GetIdent().c_str(), get_type_id(arg->GetType()) );
  }
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagInfixOprt::GenDesc() {
  string & str = this->descr;
  stdsprintf(str, "operator %s (%s)", this->GetIdent().c_str(), get_type_id( this->func_sig[2] )); 
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagInfixOprt::GetProtoType() const {
 return _T(this->func_sig.c_str());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagInfixOprt::GetDesc() const{
  return _T(this->descr.c_str());
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken* diagInfixOprt::Clone() const{
  return new diagInfixOprt(*this);
}
/*********************************************************************
  class diagOprtBin : class interface for internal binary operators
 **********************************************************************/
diagOprtBin::diagOprtBin(
    model_class * model, const char * token_str, 
    EOprtPrecedence priority, EOprtAsct assoc,
    pair<string, diag_bin_oprt_func> & flist, bool optimize
    )
     : IOprtBin( _T(token_str), (int)priority, assoc )
{
  rassert( model != NULL );
  this->mo = model;
  this->func_sig = flist.first;
  this->callback = flist.second;
  // check stack
  int argc = get_stack_argc(func_sig);
  if( optimize ){ // only allow ints, floats, complexes as arguments as types will not be checked at runtime
    assert_stack_can_be_optimized(func_sig);
  } 
  rassert( callback != NULL );
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagOprtBin::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  if(  mo->iscompiled() ){
    this->callback(ret,   (argv[0]).Get(),  (argv[1]).Get() );
    return;
  }else{
    // strict type check, etc etc ...
    assert( argc == 2 );
    const ptr_val_type & arg1 = argv[0];
    const ptr_val_type & arg2 = argv[1];
    // first we determine which function to call based on the arg stack type
    if( does_stack_match( func_sig, argv, argc ) ){
      (  callback  ) ( ret,  arg1.Get(),  arg2.Get() );
      return;
    } 
  }
  throw_prototype_error(this, argv, argc);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagOprtBin::GetProtoType() const {
 return _T(this->func_sig.c_str());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagOprtBin::GetDesc() const{
  return this->descr.c_str();
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken* diagOprtBin::Clone() const{
  return new diagOprtBin(*this);
}
/*********************************************************************
  class profileOprtBin : class interface for internal binary operators used on profiles
 **********************************************************************/
profileOprtBin::profileOprtBin(
    model_class * model, const char * token_str, 
    EOprtPrecedence priority, EOprtAsct assoc,
    pair<string, diag_bin_oprt_func> & flist, const int index 
    )
     : IOprtBin( _T(token_str), (int)priority, assoc ), pindex(index)
{
  rassert( model != NULL );
  this->mo = model;
  this->func_sig = flist.first;
  this->callback = flist.second;
  // check stack
  int argc = get_stack_argc(func_sig);
  rassert( callback != NULL );
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// Have to build some special aux functions to avoid all the complications
// of casting TokenPtr<IValue> in and out of constness ...
bool does_stack_match_special(const string & fcstack, const IValue * arg1, const IValue * arg2){
  const int clen = fcstack.size() - 2; // ignore the return and colon, ex: i:fffa -> fffa
  const char * cstack = &((fcstack.c_str())[2]);
  for(int iter = 0; iter < clen ; iter++){
    if( cstack[iter] == '*' ) break;
    if(iter > 2) return false;
    if( ! does_value_conform_to_type(
          ( iter == 0 ? *arg1 : *arg2 ), cstack[iter]) ) return false;
  }
  return true;
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void profileOprtBin::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  // strict type check, etc etc ...
  assert( argc == 2 );
  const ptr_val_type & arg1 = argv[0];
  const ptr_val_type & arg2 = argv[1];
  bool process_array_args = ( pindex >= 0 && ( ( arg1->GetType() == 'a' && arg2->GetType() != 'a' )
        ||( arg2->GetType() == 'a' && arg1->GetType() != 'a' ) ) ); 
  if(  ! process_array_args ){
    if( does_stack_match( func_sig, argv, argc ) ){
      ( callback  ) ( ret, arg1.Get(), arg2.Get() );
      return;
    }
  }else{ // one argument is an array with implicit subtyping using this->pindex value 
    const int aindex = arg1->GetType() == 'a' ? 0 : 1; // array index
    const int sindex = arg1->GetType() != 'a' ? 0 : 1; // scalar index
    const array_type & array = argv[aindex]->GetArray();
    if( this->pindex < array.size() ){
      const IValue * array_value = const_cast<Value*>(& (array[pindex]) )->AsIValue();
      // modify the func signature to match the 1st argument for the scalar if it is an array
      string modified_func_sig(this->func_sig);
      modified_func_sig[ aindex + 2 ] = modified_func_sig[ sindex + 2 ];
      if( does_stack_match_special( modified_func_sig,  array_value,  argv[sindex].Get() ) ){ 
        ( callback  ) ( ret, array_value, argv[sindex].Get() );
        return;
      }
      // array_type value and scalar arg types did NOT match any call stack arguments
      throw_prototype_error_special(this, array_value, argv[sindex].Get() ); 
    }
  }
  throw_prototype_error(this, argv, argc);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * profileOprtBin::GetProtoType() const{
  return _T(this->func_sig.c_str());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * profileOprtBin::GetDesc() const{
  return _T(this->descr.c_str());
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken* profileOprtBin::Clone() const{
  return new profileOprtBin(*this);
}
/*********************************************************************
  class builtin funcs: class that handles various overloaded builtin functions
 **********************************************************************/
diagBuiltInFunc::diagBuiltInFunc(model_class * model, const char * name, 
    pair<string, diag_builtin_func> & flist, bool optimized)
                             : ICallback( cmFUNC, name, -1 )
{
  rassert( name != NULL );
  this->mo = model;
  this->func_sig = flist.first;
  this->callback = flist.second;
  int numargs = get_stack_argc(func_sig); 
  if( optimized ) rassert( func_sig[0] == 'c'  || func_sig[0] == 'f' || func_sig[0] == 'i' );
  rassert( callback != NULL );
  this->GenDesc();
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagBuiltInFunc::GenDesc(){
  // string & str = this->descr;
  const string & name = this->GetIdent();
  string soverload;
  gen_call_desc(soverload, "", this->GetIdent(), this->func_sig);
  this->descr += soverload; 
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diagBuiltInFunc::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  if( mo->iscompiled() ){
    // just call optimizied function
    callback( ret, argv, argc ); return;
  }else{
    // find the correct overloaded function to call and if none is found throw an error
    const string & call_stack_types = this->func_sig;
    if( does_stack_match( call_stack_types, argv, argc ) ){
      ( callback ) (ret, argv, argc) ; return;
    }
  }
  throw_prototype_error(this, argv, argc, false);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagBuiltInFunc::GetProtoType() const {
  return _T(this->func_sig.c_str());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diagBuiltInFunc::GetDesc() const {
  return _T(this->descr.c_str());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken * diagBuiltInFunc::Clone() const{
  return new diagBuiltInFunc(*this);
};

/*********************************************************************
  class diag_func_class : class interface for internal diag functions
 **********************************************************************/
diag_func_class::diag_func_class(diag_func cback, const char * name, const string & call_sig,
    bool opt)
             : ICallback( cmFUNC, name, 0 )
{
  rassert(cback != NULL);
  rassert(name != NULL);
  this->optimized = opt;
  this->func_ptr = cback;
  this->name = string(name);
  this->func_sig = call_sig;
  const int argc = get_stack_argc(func_sig); // throws an error if the stack is invalid
  // usually functions are not variadic:
  this->variadic = (argc >= 0);
  this->SetArgc(argc);
  this->GenDesc();
  // now check if handles are part of the stack, if so, no optimization is allowed
  for(int iter = 2; iter < this->func_sig.size(); iter++){
    if( func_sig[iter] == 'h' ){ this->optimized = false; break; }
  }
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diag_func_class::Eval(ptr_val_type & ret, const ptr_val_type * argv, int argc){
  model_class & mo = getModel();
  bool compiled = mo.iscompiled();
  if( ( ( ! compiled ) || ( ! optimized ) ) &&  ! does_stack_match(this->func_sig, argv, argc) ){
    throw_prototype_error(this, argv, argc, false);
  }
  this->func_ptr( *this, mo, compiled, ret, argv, argc );
  if( (  ! compiled || ! optimized ) && ! does_return_match(this->func_sig, ret)){
    throw_return_error( this, ret, argv, argc);
  }
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diag_func_class::GetProtoType() const {
 return  _T(this->func_sig.c_str());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const char_type * diag_func_class::GetDesc() const {
  return _T(this->descr.c_str());
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IToken * diag_func_class::Clone() const{
  return new diag_func_class(*this);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void diag_func_class::GenDesc(){
  gen_call_desc(this->descr, "", this->name, this->func_sig);
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
const string & diag_func_class::GetName() const{
  return this->name;
};
/*   
     end class diag_func_class 

 ********************************************************************************/

//-----------------------------------------------------------------------------------//
// forward declaration
int diagprintf(const qstream &, const char *, const diag_func_class &, const ptr_val_type *, const int);
//-----------------------------------------------------------------------------------//

void resetdiagout(){
  if(!!stored_diagout.is_closed()){ 
    stored_diagout.close();
  }
};

void setdiagout(const qstream & ndiag){
  stored_diagout = ndiag;
};

const qstream & getdiagout(){ 
  if( ! stored_diagout.is_closed() ) return stored_diagout;
  return getModel().getSTDOUT();
};
//-----------------------------------------------------------------------------------//
#define DIAGFUNC(name) void name(diag_call_stack)
#define getarg(iarg) (*(argv[(iarg)].Get()))
#define defarg(token,iarg) const Value & token = (*(argv[(iarg)].Get())); rassert(iarg < argc);
#define RETURNDIAG(expr) { ( *ret = (expr) ); return; };
#define RETURNIFNOTCOMPILED(val) { if(! compiled ){ RETURNDIAG(val); } }
#define __thisdiagfunc__ ( func_class.GetName().c_str() )
//-----------------------------------------------------------------------------------//
/* Global Array Indexing function : returns a pointer to array[index]  */
Variable * diagOprtArrayIndex( const array_type & Array, const IValue & Index ){
  static Value dummy(0); // return this if compilation is not done 
  if(! does_value_conform_to_type( Index, 'i' ) ){ // bad index type
    mpError("Array indexing operator [] must take an integer as its argument.");
  }
  if(  getModel().iscompiled() ){
    int_type i = Index.GetInteger();
    if( i >= 0 && i < Array.size() ){
      return new Variable( const_cast<Value*>(&Array[i]) );
    }else{
      mpError("Array accessed out of bounds.");
    }
  }
  return new Variable(& dummy );
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_rewindtime) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  int_type numsteps = getarg(0).GetInteger();
  psystem & p = mo.getPsystem();
  if(numsteps <= 0 || numsteps >= p.get_tindex() + 1){
    RETURNDIAG(rval=false);
  }
  unsigned int nsteps = static_cast<unsigned int>(numsteps);
  p.rewindtime(nsteps);
  RETURNDIAG(rval=true);
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_settimestep) {
  bool rval;
  RETURNIFNOTCOMPILED(rval=false);
  float_type new_tstep = getarg(0).GetFloat();
  REAL tstep = static_cast<REAL>(new_tstep);
  if(tstep <= 0) RETURNDIAG(rval=false);
  int iter = mo.get_global_param_index_from_name(gdb_str_TIMESTEP);
  parameter & tstepparam = mo.get_global_param_by_index(iter);
  tstepparam.setValue(tstep);
  RETURNDIAG(rval=true);
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_getprevtimestep) {
  float_type rvalue;
  RETURNIFNOTCOMPILED(rvalue=0.0);
  int_type tindex = getarg(0).GetInteger();
  if(tindex < 0) mpError("tindex argument to %s must be greater than or equal to 0.", __thisdiagfunc__);
  ssystem & s = mo.getSsystem();
  psystem & p = mo.getPsystem();
  unsigned int itindex = static_cast<unsigned int>(tindex);
  if(itindex > static_cast<unsigned int>(p.get_tindex())){
    mpError("tindex argument to '%s' exceeds the current number of time steps.", __thisdiagfunc__);
  }
  RETURNDIAG(rvalue=s.gettimestep(itindex));
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_gettimestep){
  float_type rvalue;
  REAL tstep = mo.get_global_param_by_name(gdb_str_TIMESTEP).getfValue();
  RETURNDIAG( rvalue = tstep );
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_gettime) {
  float_type rvalue = 0;
  RETURNIFNOTCOMPILED(rvalue = 0);
  int_type tindex = getarg(0).GetInteger();
  if(tindex < 0) mpError("tindex argument to '%s' must be greater than or equal to 0.", __thisdiagfunc__);
  ssystem & s = mo.getSsystem();
  psystem & p = mo.getPsystem();
  unsigned int itindex = static_cast<unsigned int>(tindex);
  if(itindex > static_cast<unsigned int>(p.get_tindex())){
    mpError("tindex argument to '%s' exceeds the current number of time steps.", __thisdiagfunc__);
  }
  RETURNDIAG(rvalue=s.getelapsedtime(itindex));
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_getnumtsteps){
  int_type rvalue;
  RETURNIFNOTCOMPILED(rvalue=-1);
  psystem & p = mo.getPsystem();
  int nsteps = p.getnumtsteps();
  RETURNDIAG( rvalue=nsteps );
}
//-----------------------------------------------------------------------------------//
DIAGFUNC(diag_func_tindex){
  int_type rvalue;
  RETURNIFNOTCOMPILED(rvalue=-1);
  psystem & p = mo.getPsystem();
  RETURNDIAG( rvalue = p.get_tindex() );
}
// ............................................................... //
enum type_S_data { phi = 1, G = 2, R = 3, S = 4, Fphi = 5, FG = 6, none };
// ............................................................... //
type_S_data convertint2type_S_data(const int & type, const diag_func_class & func_class, 
    bool throwexcepinvalid = true){
  // use good c++ here,  no direct casts!
  switch(static_cast<int>(type)){
    case 1 : return(phi);
    case 2 : return(G);
    case 3 : return(R);
    case 4 : return(S);
    case 5 : return(Fphi);
    case 6 : return(FG);
    default:   
             if(throwexcepinvalid){
               mpError("Invalid S-system data set provided to function '%s'.", __thisdiagfunc__);
             }
             return none;
  };
};
// ............................................................... //
const regset & diag_get_regset_aux(const type_S_data & datatype, ssystem & ssys, int t, int tprime, bool
    return_regularized_set){
  switch(datatype){
    case phi:
      return(ssys.get_phi_set(t, tprime));
      break;
    case G:
      return(ssys.get_G_set(t, tprime));
      break;
    case R:
      if(return_regularized_set) return(ssys.get_R_regularized_set(t));
      return(ssys.get_R_set(t, tprime));
      break;
    case S:
      if(return_regularized_set) return(ssys.get_S_regularized_set(t));
      return(ssys.get_S_set(t, tprime));
      break;
    case Fphi:
      if(return_regularized_set) return(ssys.get_F_phi_regularized_set(t));
      return(ssys.get_F_phi_set(t, tprime));
      break;
    case FG:
      return(ssys.get_F_G_set(t, tprime));
      break;
    default:
      fatal_stop("Invalid register type requested.");
      return(ssys.get_F_G_set(t, tprime)); // avoids a warning
  };
};
// ............................................................... //
const qgcomplex diag_func_get_reg_data_aux(
    const type_S_data & datatype, 
    bool provisional,
    const diag_func_class & func_class,
    const ptr_val_type * argv, const int argc)
{
  qgcomplex zero(0,0);
  model_class & mo = getModel();
  ssystem & ssys = mo.getSsystem();
  psystem & p = mo.getPsystem();
  int t, tprime, ki, nki, i, n;
  if(argc > 5){
    t      = static_cast<int>( getarg(0).GetInteger() );
    tprime = static_cast<int>( getarg(1).GetInteger() );
    nki    = static_cast<int>( getarg(2).GetInteger() );
    ki     = static_cast<int>( getarg(3).GetInteger() );
    i      = static_cast<int>( getarg(4).GetInteger() );
    n      = static_cast<int>( getarg(5).GetInteger() );
  }else{
    tprime = static_cast<int>( getarg(0).GetInteger() );
    t      = tprime;
    nki    = static_cast<int>( getarg(1).GetInteger() );
    ki     = static_cast<int>( getarg(2).GetInteger() );
    i      = static_cast<int>( getarg(3).GetInteger() );
    n      = static_cast<int>( getarg(4).GetInteger() );
  }
  const int ksize = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData().size();
  const int nksize = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData().size();
  const int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  if(p.get_tindex() >= 0){
    if(t < tprime){
      mpError("Time index argument t' must be less than or equal to t in '%s'.", __thisdiagfunc__);
    }
    if(tprime < 0 || t < 0){
      mpError("Time index arguments must be greater than or equal to 0 in '%s'.", __thisdiagfunc__);
    }
    if(tprime > p.get_tindex() || tprime > p.get_tindex()){
      mpError("Time index arguments must be less than or equal to %i in '%s'.",p.get_tindex(), __thisdiagfunc__);
    }
    if(nki >= nksize || nki <= -nksize){
      mpError( "Fourier angular index must lie between '%i' and '%i' in '%s'.", 0, nksize-1, __thisdiagfunc__);
    }
    if(ki < 1 || ki > ksize){
      mpError( "Fourier mod 'ki' must lie between '%i' and '%i' in '%s'.",1, ksize,__thisdiagfunc__);
    }
    if(i < 1 || i >numfields || n < 1 || n > numfields){
      mpError( "Field index arguments must lie between '%i' and '%i' in '%s'",1,numfields,__thisdiagfunc__);
    }
    regset reg = diag_get_regset_aux(datatype, ssys, t, tprime, false);
    if(provisional){
      reg.setintermediate();
    }else{
      reg.setprimary();
    }
    const qgcomplex & val = ssys.get_register( reg, nki, ki-1, i-1, n-1 );
    // if we have requested an equal time set which has a regularized component and the model is using regularization
    // then add to this value the regularized piece
    if( t==tprime && ssys.is_model_using_regularization() && ( datatype == R || datatype == S || datatype == Fphi )){
      regset regularized_reg = diag_get_regset_aux(datatype, ssys, t, tprime, true);
      if(provisional){
        regularized_reg.setintermediate();
      }else{
        regularized_reg.setprimary();
      };
      const qgcomplex regularized_val = ssys.get_register( regularized_reg, nki, ki-1, i-1, n-1 );
      return( regularized_val + val );
    }
    return(val);
  }
  return zero;
};
// ............................................................... //
DIAGFUNC(diag_func_error) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  // basically just call printf function and throw an error
  const string & format = getarg(0).GetString();
  diagprintf(getdiagout(), format.c_str(), func_class, &(argv[1]), argc-1); 
  mpError("\nUser generated error has requested program termination.");
  RETURNDIAG(rval=false);
}
// ............................................................... //
int getrequireddigits(double val){
  // get the integer part first
  int ndigits = 1;
  if(val < 0){
    val = -val;
    ndigits++;
  };
  double ival = static_cast<double>(static_cast<int>(val)); // truncate
  val=val-ival;
  if(ival > 0){
    ndigits += static_cast<int>(log10(ival));
  }else{
    ndigits++;
  }
  while( ndigits < 40 ){ 
    if( val-static_cast<double>(static_cast<int>(val)) == 0.0) break;
    val*=10.0;
    ndigits++;
  }
  return(ndigits);
};

// ............................................................... //
DIAGFUNC(diag_func_getmaxdigits) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  int dsize = 0;
  for(int iter = 0; iter < argc; iter++){
    int rsize = getrequireddigits(getarg(0).GetFloat());
    if(rsize > dsize) dsize = rsize;
  }
  RETURNDIAG(rval=dsize);
}
// ............................................................... //
DIAGFUNC(diag_func_bessel1) {
  float_type rval;
  rval = static_cast<float_type>( j1( getarg(0).GetFloat() ) );
  RETURNDIAG(rval);
}
// ............................................................... //
DIAGFUNC(diag_func_set_static_data) {
  bool_type rval;
  const int dataset = static_cast<int>( getarg(0).GetInteger() );
  if(dataset != 0 && dataset != 1){
    mpError("Invalid static data specified as first argument of function %s.", __thisdiagfunc__);
  }
  RETURNIFNOTCOMPILED(rval=false);
  //const char * fmt = "(dataset, complex value, nki, npi, k, p, q, i, j, k, m, a, b)";
  ssystem & ssys = mo.getSsystem();
  static_container * coeffs = dataset == 0 ? ssys.getRcoeffs_dyn() : ssys.getScoeffs_dyn();
  const cmplx_type & setvalue = getarg(1).GetComplex();
  rank11 r11;
  r11[0] = static_cast<INDEX>(getarg(2).GetInteger());
  r11[1] = static_cast<INDEX>(getarg(3).GetInteger());
  r11[2] = static_cast<INDEX>(getarg(4).GetInteger()-1);
  r11[3] = static_cast<INDEX>(getarg(5).GetInteger()-1);
  r11[4] = static_cast<INDEX>(getarg(6).GetInteger()-1);
  r11[5] = static_cast<INDEX>(getarg(7).GetInteger()-1);
  r11[6] = static_cast<INDEX>(getarg(8).GetInteger()-1);
  r11[7] = static_cast<INDEX>(getarg(9).GetInteger()-1);
  r11[8] = static_cast<INDEX>(getarg(10).GetInteger()-1);
  r11[9] = static_cast<INDEX>(getarg(11).GetInteger()-1);
  r11[10]= static_cast<INDEX>(getarg(12).GetInteger()-1);
  try{
    coeffs->storevaluediag(r11, qgcomplex(setvalue.real(), setvalue.imag()));
    RETURNDIAG(rval=true);
  }catch(const cgcmexception & e){
    RETURNDIAG(rval=false);
  }
};
// ............................................................... //
DIAGFUNC(diag_func_get_static_data) {
  const int dataset = static_cast<int>( getarg(0).GetInteger() );
  if(dataset != 0 && dataset != 1){
    mpError("Invalid static data specified as first argument of function %s.", __thisdiagfunc__);
  }
  RETURNIFNOTCOMPILED(cmplx_type(1,1));
  //const char * fmt = "(dataset, nki, npi, k, p, q, i, j, k, m, a, b)";
  ssystem & ssys = mo.getSsystem();
  const static_container * coeffs = dataset == 0 ? ssys.getRcoeffs_dyn() : ssys.getScoeffs_dyn();
  rank11 r11;
  r11[0] = static_cast<INDEX>(getarg(1).GetInteger());
  r11[1] = static_cast<INDEX>(getarg(2).GetInteger());
  r11[2] = static_cast<INDEX>(getarg(3).GetInteger()-1);
  r11[3] = static_cast<INDEX>(getarg(4).GetInteger()-1);
  r11[4] = static_cast<INDEX>(getarg(5).GetInteger()-1);
  r11[5] = static_cast<INDEX>(getarg(6).GetInteger()-1);
  r11[6] = static_cast<INDEX>(getarg(7).GetInteger()-1);
  r11[7] = static_cast<INDEX>(getarg(8).GetInteger()-1);
  r11[8] = static_cast<INDEX>(getarg(9).GetInteger()-1);
  r11[9] = static_cast<INDEX>(getarg(10).GetInteger()-1);
  r11[10]= static_cast<INDEX>(getarg(11).GetInteger()-1);
  try{
    const qgcomplex & val = coeffs->getvalue(r11);
    RETURNDIAG( cmplx_type( val.real(), val.imag() )  );
  }catch(const cgcmexception & e){
    RETURNDIAG( cmplx_type(0,0) );
  }
};
// ............................................................... //
DIAGFUNC(diag_func_string) {
  // takes variable converts to string if possible
  string_type rval;
  RETURNIFNOTCOMPILED(rval);
  rval = cast_to_string( getarg(0) );
  RETURNDIAG( rval );
};
// ............................................................... //
DIAGFUNC(diag_func_real) {
  // takes variable converts to real if possible
  float_type rval;
  RETURNIFNOTCOMPILED(rval=0.1);
  rval = cast_to_float( getarg(0) );
  RETURNDIAG( rval );
};
// ............................................................... //
DIAGFUNC(diag_func_imag) {
  // takes variable converts to real if possible
  float_type rval;
  const cmplx_type & cmplx = getarg(0).GetComplex();
  RETURNDIAG( rval=(cmplx.imag()) );
};
// ............................................................... //
DIAGFUNC(diag_func_integer) {
  // takes variable converts to real if possible
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  rval = cast_to_int( getarg(0) );
  RETURNDIAG( rval );
};
// ............................................................... //
DIAGFUNC(diag_func_bool) {
  // takes variable converts to real if possible
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  rval = cast_to_bool( getarg(0) );
  RETURNDIAG( rval );
};
// ............................................................... //
DIAGFUNC(diag_func_complex) {
  // takes variable converts to real if possible
  cmplx_type rval(1,1);
  RETURNIFNOTCOMPILED(rval);
  rval = cmplx_type( getarg(0).GetFloat(), getarg(1).GetFloat() );
  RETURNDIAG( rval );
};
// ............................................................... //
DIAGFUNC(diag_func_getSsysdata) {
  // check number of args
  if(argc < 7){
    mpError("Function %s takes at least 7 arguments: %s(data set id, bool provisional, t, [t'], nki, ki, i,,n).",
        __thisdiagfunc__,__thisdiagfunc__);
  }
  RETURNIFNOTCOMPILED(cmplx_type(1,1));
  type_S_data type = convertint2type_S_data(getarg(0).GetInteger(), func_class, true); // fail if not valid type
  bool provisional = getarg(1).GetBool();
  const ptr_val_type * subargv = &(argv[2]);
  const int subargc = argc - 2;
  const qgcomplex cval = diag_func_get_reg_data_aux(type, provisional, func_class, subargv, subargc);
  RETURNDIAG(cmplx_type(cval.real(), cval.imag()));
}
// ............................................................... //
const char * getEscapeSequence(char c){
  char sequences[12] = { '\'', '\"', '\\', '0', 'a', 'b', 'f', 'n', 'r', 't',
    'v', '%' };
  const char * escsequences[12] = {
    "\'", "\"", "\\", "\0", "\a", "\b", "\f", "\n", "\r", "\t", "\v", "%" };
  for(int i = 0; i < 12; i++){
    if(c == sequences[i]) return(escsequences[i]);
  }
  return("\0");
}
// ............................................................... //
DIAGFUNC(diag_func_fclose) {
  int_type rval = 0;
  for(int i = 0; i < argc; i++){
    // check types, make sure they are handles
    if( is_value_valid_handle( getarg(i) ) ){
      if( compiled ){
        if(handler.file_close( getarg(i) )) rval++;
      }
    }
  }
  RETURNDIAG(rval);
};
// ............................................................... //
DIAGFUNC(diag_func_fopen) {
  RETURNIFNOTCOMPILED(get_default_value_of_type('h'));
  const string & fname = getarg(0).GetString();
  const string & opt = getarg(1).GetString();
  RETURNDIAG( handler.file_open( fname.c_str(), opt.c_str() ));
};
// ............................................................... //
DIAGFUNC(diag_func_getstdout) {
  RETURNIFNOTCOMPILED(get_default_value_of_type('h'));
  RETURNDIAG( handler.getstdout()  );
};
// ............................................................... //
DIAGFUNC(diag_func_createarray){
  // fully variadic
  RETURNIFNOTCOMPILED(get_default_value_of_type('a'));
  Value array( argc, static_cast<int_type>(0) );
  for(int i = 0; i < argc; i++){
    array[i] = getarg(i);
  }
  RETURNDIAG(array);
};
// ............................................................... //
DIAGFUNC(diag_func_arrayelement){
  RETURNIFNOTCOMPILED(get_default_value_of_type('v'));
  const array_type & array = getarg(0).GetArray();
  int_type index = getarg(1).GetInteger();
  if( index >= 0 && index < array.size() ){
    RETURNDIAG( array[index] );
  }
  mpError("second argument of function '%s' must be valid array index for first argument.", __thisdiagfunc__);
  RETURNDIAG(0);
};
// ............................................................... //
bool ischarmatch(char c, const char * str){
  if(str){
    for(size_t i = 0; i < strlen(str); i++){
      if(str[i] == c) return(true);
    }
  }
  return( false);
};
// ............................................................... //
int diagprintfvar(const qstream & fout, const IValue & pvalue, char fmtspec, const string & flags,
    int width, int precision){
  string fmt = "%" + flags;
  if(width > 0) fmt+= convertint2string(width);
  if(precision > 0) fmt+= string(".") + convertint2string(precision);
  fmt+=string(" ");
  fmt[fmt.size()-1] = fmtspec;
  if(fmtspec == 'c'){
    return fout.printf(fmt.c_str(), cast_to_char(pvalue));
  }else if(ischarmatch(fmtspec, "diouxXp")){
    int ivalue = cast_to_int(pvalue);
    return fout.printf( fmt.c_str(),ivalue);
  }else if(ischarmatch(fmtspec, "eEfgG")){
    float fvalue = cast_to_float(pvalue);
    return fout.printf( fmt.c_str(), fvalue);
  }else if(fmtspec == 's'){
    // a string , requires conversion
    string data = cast_to_string(pvalue);
    return fout.printf( fmt.c_str(), data.c_str());
  }else if(fmtspec == '%'){
    return fout.printf( "%%");
  }
  return 0;
};
// ............................................................... //
int readdiagprintfflags(string & sformat, const char * format, int spos){
#ifdef DEBUG
  assert(spos < static_cast<int>(strlen(format)));
#endif
  if( ischarmatch( format[spos], "-+ #") ){
    sformat += string(" ");
    sformat[ sformat.size() - 1] = format[spos];
    return( 1 );
  }
  return( 0);
};
// ............................................................... //
int readdiagprintfwidth(const diag_func_class & func_class, const char * format, const int spos, 
    const ptr_val_type * argv, const int list_size, int & arg_pos, int & width){
  int pos = spos;
  if( format[pos] == '*' ){
    if(arg_pos >= list_size){ 
      mpError("Function '%s' expects more arguments than were supplied.", __thisdiagfunc__); 
    };  
    width = cast_to_int( getarg(arg_pos) ); //static_cast<int>( list[arg_pos] );
    pos++; arg_pos++;
    if(width <= 0){
      mpError("Function '%s' expects positive width/precision argument.", __thisdiagfunc__);
    }
  }else if(ischarmatch(format[pos], "0123456789")){
    int shift = 1;
    while( format[pos] >= '0' && format[pos] <= '9'){
      if(width < 0) width = 0;
      width *= shift;
      width += static_cast<int>(( format[pos] - '0' ));
      shift *= 10;
      pos++;
    };
    if(width <= 0){
      mpError("Function '%s' expects positive width/precision argument.", __thisdiagfunc__);
    }
  }
  return( pos-spos);
};
// ............................................................... //
int readdiagprintfprec(const diag_func_class & func_class, const char * format, const int spos,
    const ptr_val_type * list, const int list_size, int & arg_pos, int & precision){
  if( format[spos] != '.' ) return( 0);
  return( readdiagprintfwidth(func_class, format, spos+1, list, list_size, arg_pos, precision)+1);
};
// ............................................................... //
int readdiagprintflengthspec(const char * format, const int spos){
  //if(format[spos] == 'h' || format[spos] == 'l' || format[spos] == 'L'){
  if(ischarmatch( format[spos], "hlL" )){
    return( 1);
  }
  // this function simply reads these specifiers, they are not actually used
  return( 0);
};
// ............................................................... //
int readfmtspec(const diag_func_class & func_class, const char * format, const int spos, char & fmtspec)
{
  int pos = spos; 
  if( ischarmatch( format[pos], "cdieEfgGosuxXp%" ) ){
    fmtspec = format[pos];
  }else{
    mpError("Unrecognized format specifier '%c' supplied in function '%s'.", format[pos],__thisdiagfunc__);
  }
  return( pos - spos);
}; 
// ............................................................... //
int diagprintf(const qstream & fout, const char * format, const diag_func_class & func_class, 
    const ptr_val_type * argv, const int argc)
{
  int rvalue = 0;
  int arg_pos = 0;
  int fmt_pos = 0;
  char cf = format[0];
  while( cf != 0 ){
    if( cf == '%' || cf == '\\' ){ // we have a format identifier or escape sequence
      // flush the stream
      if( cf == '\\'){
        fmt_pos++; if(format[fmt_pos] == 0) break; // end of format
        const char * sequence = getEscapeSequence(format[fmt_pos]);
        rvalue+=fout.printf("%s", sequence);
      }else{ // '%' it's a format specifier
        int precision = -1; int width = -1; 
        string flags; 
        char fmtspec;
        fmt_pos++; if(format[fmt_pos] == 0) break;
        fmt_pos+=readdiagprintfflags(flags, format, fmt_pos); 
        if(format[fmt_pos] == 0) break;
        fmt_pos+=readdiagprintfwidth(func_class, format, fmt_pos, argv, argc, arg_pos, width); 
        if(format[fmt_pos] == 0) break;
        fmt_pos+=readdiagprintfprec (func_class, format, fmt_pos, argv, argc, arg_pos, precision); 
        if(format[fmt_pos] == 0) break;
        fmt_pos+=readdiagprintflengthspec(format, fmt_pos); 
        if(format[fmt_pos] == 0) break;
        fmt_pos+=readfmtspec(func_class,format, fmt_pos, fmtspec); 
        if(arg_pos >= argc){
          mpError("Function %s expects more arguments than were supplied.", __thisdiagfunc__); 
        }
        const IValue & pvalue = getarg(arg_pos); arg_pos++;
        rvalue += diagprintfvar(fout, pvalue, fmtspec, flags, width, precision);
        if(format[fmt_pos] == 0) break; 
      }
    }else{ 
      // write character to the stream
      if(fout.putc(cf)) rvalue++;
      //if(fputc( static_cast<int>( cf ), fout ) == cf) rvalue++;
    };
    fmt_pos++; cf = format[fmt_pos];
  }
  return rvalue;
};
// ............................................................... //
DIAGFUNC(diag_func_printf) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  // get the fmt string
  const string & fmt = getarg(0).GetString();
  rval = diagprintf(getdiagout(), fmt.c_str(), func_class, &(argv[1]), argc-1);
  RETURNDIAG(rval); 
}
// ............................................................... //
DIAGFUNC(diag_func_flush) { // flush the file handle
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  for(int iter = 0; iter < argc; iter++){
    if( ! is_value_valid_handle(getarg(iter) ) ) RETURNDIAG( rval=false );
    const qstream * fout = handler.get_file( getarg(0) );
    if(fout != NULL){
      fout->flush();
    }
  }
  RETURNDIAG(rval=true);
}
// ............................................................... //
DIAGFUNC(diag_func_print) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  const qstream & fout = getdiagout();
  for(int iter = 0; iter < argc; iter++){
    const IValue & arg = getarg(iter);
    fout << arg;
    rval++;
  }
  fout.printf("\n");
  RETURNDIAG(rval);
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getA) {
  float_type rval;
  RETURNIFNOTCOMPILED(rval=0.0);
  // returns the actual (p x q)_3 * A_{ijk}(k,p,q)
  const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
  // first 3 args are ki, pi, qi
  INDEX ki = static_cast<INDEX>( getarg(0).GetInteger() - 1 );
  INDEX pi = static_cast<INDEX>( getarg(1).GetInteger() - 1 );
  INDEX qi = static_cast<INDEX>( getarg(2).GetInteger() - 1 );
  if( ki < 0 || ki >= static_cast<INDEX>(kc.size()) 
      || pi < 0 || pi >= static_cast<INDEX>(kc.size()) 
      || qi < 0 || qi >= static_cast<INDEX>(kc.size())){
    mpError("Kindex arguments to function '%s' must lie within range of defined domain.", __thisdiagfunc__);
  }
  const int nfs = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  INDEX i = static_cast<INDEX>( getarg(3).GetInteger() );
  INDEX j = static_cast<INDEX>( getarg(4).GetInteger() );
  INDEX k = static_cast<INDEX>( getarg(5).GetInteger() );
  if( i < 1 || i > nfs || j < 1 || j > nfs || k < 1 || k > nfs ){
    mpError("Field index arguments must lie between 1 and %i.\nBad field index argument in function '%s'.",
        nfs, __thisdiagfunc__);
  }
  // now the field index pair must exist within the field index set, or we return 0
  // this function should return \hat{A} not the actual A operator, for this reason, we introduce
  // we have changed it to avoid returning infinite values where inappropriate 
  REAL rvalue = mo.computeA(rank6(ki, pi, qi, i - 1,j- 1,k - 1) );
  RETURNDIAG( rval = rvalue );
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getk) {
  float_type rval;
  RETURNIFNOTCOMPILED(rval=0.0);
  const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
  INDEX ki = static_cast<int>(getarg(0).GetInteger());
  if( ki < 1 || ki > static_cast<INDEX>(kc.size()) ){
    mpError("Kindex arguments to function '%s' must lie within range of defined domain.", __thisdiagfunc__);
  }
  REAL fval = kc[ki-1];
  RETURNDIAG(rval=fval);
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getnk) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  const vector<int> & nkdom = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData();
  INDEX nki = static_cast<int>(getarg(0).GetInteger());
  if( nki < 0 || nki >= static_cast<INDEX>(nkdom.size()) ){
    mpError("nkindex arguments to function '%s' must lie within range of defined domain.", __thisdiagfunc__);
  }
  int ival = nkdom[nki];
  RETURNDIAG(rval=ival);
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getnumks) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0.0);
  const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
  RETURNDIAG(rval = kc.size() )
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getdk) {
  float_type rval;
  RETURNIFNOTCOMPILED(rval=0.0);
  const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_DK_DOMAIN).getData();
  INDEX ki = static_cast<int>(getarg(0).GetInteger());
  if( ki < 1 || ki > static_cast<INDEX>(kc.size()) ){
    mpError("Kindex arguments to function '%s' must lie within range of defined domain.", __thisdiagfunc__);
  }
  REAL fval = kc[ki-1];
  RETURNDIAG(rval=fval);
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getQ) {
  float_type rval;
  RETURNIFNOTCOMPILED(rval=0.0);
  const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
  INDEX i = static_cast<INDEX>(getarg(1).GetInteger());
  INDEX j = static_cast<int>(getarg(2).GetInteger());
  INDEX ki = static_cast<int>(getarg(0).GetInteger() - 1);
  const int nfs = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  if(i < 1 || i > nfs || j < 1 || j > nfs){
    mpError("Field index arguments must lie between 1 and %i.\nBad field index argument in function '%s'.",
        nfs, __thisdiagfunc__);
  }
  if( ki < 0 || ki >= static_cast<INDEX>(kc.size()) ){
    mpError("Kindex arguments to function '%s' must lie within range of defined domain.", __thisdiagfunc__);
  }
  REAL fval = mo.computeQ(rank3(ki,i-1, j-1) );
  RETURNDIAG(rval=fval);
};
// ------------------------------------------------------------------------------------------- //
DIAGFUNC(diag_func_getQI) {
  float_type rval;
  RETURNIFNOTCOMPILED(rval=0.0);
  const vector<REAL> & kc = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData();
  INDEX i = static_cast<INDEX>(getarg(1).GetInteger());
  INDEX j = static_cast<int>(getarg(2).GetInteger());
  INDEX ki = static_cast<int>(getarg(0).GetInteger() - 1);
  const int nfs = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  if(i < 1 || i > nfs || j < 1 || j > nfs){
    mpError("Field index arguments must lie between 1 and %i.\nBad field index argument in function '%s'.",
        nfs, __thisdiagfunc__);
  }
  if( ki < 0 || ki >= static_cast<INDEX>(kc.size()) ){
    mpError("Kindex arguments to function '%s' must lie within range of defined domain.", __thisdiagfunc__);
  }
  REAL fval = mo.computeQI(rank3(ki,i-1, j-1) );
  RETURNDIAG(rval=fval);
};
//.............................................................. //
DIAGFUNC(diag_func_fprintf) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  const qstream * fout = handler.get_file(getarg(0));
  const string & fmt = getarg(1).GetString();
  if(fout != NULL){
    rval=diagprintf(*fout, fmt.c_str(), func_class, &(argv[2]), argc-2);
  }
  RETURNDIAG(rval);
}
// ............................................................... //
DIAGFUNC(diag_func_resetdiagout){
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  resetdiagout();
  RETURNDIAG(rval=true);
};
// ............................................................... //
DIAGFUNC(diag_func_system) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0);
  string cmd = getarg(0).GetString();
  rval = system(cmd.c_str());
  RETURNDIAG(rval);
}
// ............................................................... //
DIAGFUNC(diag_func_close_debug_io){
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  mo.closeDBGIO();
  RETURNDIAG(rval=true);
}
// ............................................................... //
DIAGFUNC(diag_func_set_debug_io) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  string fname = getarg(0).GetString();
  rval = mo.setDBGOUT(fname.c_str());
  rval = mo.setDBGIN(fname.c_str());
  RETURNDIAG(rval);
};
// ............................................................... //
DIAGFUNC(diag_func_setmodelstdout) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  string fname = getarg(0).GetString();
  rval = mo.setSTDOUT(fname.c_str());
  RETURNDIAG(rval);
};
// ............................................................... //
DIAGFUNC(diag_func_setmodelstdmsg) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  const string & fname = getarg(0).GetString();
  rval = mo.setSTDMSG(fname.c_str());
  RETURNDIAG(rval);
};
// ............................................................... //
DIAGFUNC(diag_func_setdiagout) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  resetdiagout();
  bool append = false;
  bool ifnotexist = false;
  // only option allowed really for now haha
  const string & fname = getarg(0).GetString();
  for(int i = 1; i < argc; i++){
    const IValue & arg = getarg(i);
    if( ! does_value_conform_to_type(arg, 's')){
      const string & opt = arg.GetString();
      if(opt == "append" || opt == "a" || opt == "APPEND") append = true;
      if(opt == "ifnotexist" || opt == "new" || opt == "n") ifnotexist = true;
    }else{
      mpError("Function %s expects argument %i to be type string.", __thisdiagfunc__, i+1, get_type_id('s'));
    }
  }
  if(ifnotexist){ // only open file if it does not exist already
    qstream newdiag;
    if(newdiag.open(fname.c_str(), "r")){
      newdiag.close(); 
      RETURNDIAG(rval=false);
    }
  }
  string opts("w");
  if(append) opts="a";
  qstream newdiag;
  rval = newdiag.open(fname.c_str(), opts.c_str());
  if(rval){
    setdiagout(newdiag);
  }
  RETURNDIAG(rval);
};
// ............................................................... //
DIAGFUNC(diag_func_ncopen) {
  int_type rval;
  RETURNIFNOTCOMPILED(rval=0); // 0 is a valid handle (i.e. NULL)
  // format:
  // filename
  // options: 0 for readonly (useless), 1 for write (newfile), 2 for overwrite
  NcFile::FileMode fmode = NcFile::ReadOnly;
  const string & filename = getarg(0).GetString();
  int_type options = getarg(1).GetInteger();
  if(options == 0){ fmode = NcFile::ReadOnly; }
  else if(options == 1){ fmode = NcFile::Write; }
  else if(options == 2){ fmode = NcFile::Replace; }
  else{
    mpError("Error in function '%s'. Invalid netcdf file mode specified.",__thisdiagfunc__);
    RETURNDIAG(rval=0);
  }
  RETURNDIAG( handler.ncfile_open( filename.c_str(), fmode ) );
};
// ............................................................... //
DIAGFUNC(diag_func_ncsetmissingvalue) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  REAL fval = static_cast<REAL>( getarg(0).GetFloat() );
  RETURNDIAG(fval=getsetmissing_val(fval, true));
}
// ............................................................... //
DIAGFUNC(diag_func_ncwriteprofiles) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  NcFile * ncfile = handler.get_nc_file(getarg(0));
  if(!ncfile) RETURNDIAG(rval=false);
  // iterate through profiles and write them to the file
  NcDim * freeindex = ncfile->get_dim("time");
  if(! freeindex->is_valid() ){
    set_mp_errmsg("Error in function '%s'. Missing necessary dimensions in file.",__thisdiagfunc__);
    RETURNDIAG(rval=false);
  }
  // seems like some definitions below are not needed
  NcDim * levelindex = ncfile->get_dim("level");
  NcDim * kindex = ncfile->get_dim("kindex");
  NcDim * nkindex= ncfile->get_dim("nkindex"); 
  int numprofs = mo.get_num_global_profs();
  for(int x = 0; x < numprofs; x++){
    const parameter & cprof = mo.get_global_prof_by_index(x);
    NcDim * pdim = NULL;
    if( (strcmp(cprof.GetName(), gdb_prof_str_K_DOMAIN_CENTERED) == 0)
        || (strcmp(cprof.GetName(), gdb_prof_str_K_DOMAIN_CENTERED_SQUARED) == 0)){ 
      pdim = kindex;
    }else if(strcmp(cprof.GetName(), gdb_prof_str_NK_DOMAIN)==0){
      pdim = nkindex;
    }
    if(cprof.getType() == parameter::String){
      // do nothing 
    }else if(cprof.getType() == parameter::Integer){
      const vector<int> & idata = cprof.getiData();
      if(static_cast<int>(idata.size()) == levelindex->size()){
        pdim = levelindex;
      }
      if(pdim==NULL){
        // add dimension with variable name
        pdim = ncfile->add_dim(cprof.GetName(), idata.size());
      }
      NcVar * var = ncfile->add_var(cprof.GetName(), ncInt, pdim); 
      var->add_att("missing_value", static_cast<int>(getsetmissing_val(0,false)));
      int * iarray = new int[idata.size()];
      for(unsigned int i = 0; i < idata.size(); i++){
        iarray[i] = idata[i];
      }
      var->put(iarray, idata.size(), 0, 0, 0, 0); 
      delete [] iarray;
    }else{
      const vector<REAL> & data = cprof.getData();
      if(static_cast<int>(data.size()) == levelindex->size()){
        pdim = levelindex;
      }
      if(pdim==NULL){
        // add dimension with variable name
        pdim = ncfile->add_dim(cprof.GetName(), data.size());
      }
      NcVar * var = ncfile->add_var(cprof.GetName(), ncREAL, pdim); 
      var->add_att("missing_value", getsetmissing_val(0,false));
      REAL * farray = new REAL[data.size()];
      for(unsigned int i = 0; i < data.size(); i++){
        farray[i] = data[i];
      }
      var->put(farray, data.size(), 0, 0, 0, 0); 
      delete [] farray;
    }
  } 
  RETURNDIAG(rval=true);
}
// ............................................................... //
DIAGFUNC(diag_func_save_model_state) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  const string & fname = getarg(0).GetString();
  const string & opts = getarg(1).GetString();
  bool overwrite = false;
  if(opts == "overwrite" || opts == "o"){
    overwrite = true;
  }
  if(overwrite){
    string aname(fname);
    // aname+=".nc"; (let the user decide)
    int rstatus = remove( aname.c_str() );
    if(rstatus != 0){ // remove failed
      // if failed because file did not exist then we are OK, otherwise return false;
      if( errno != ENOENT ){
        RETURNDIAG(rval=false);
      }
    }
  } 
  bool status = mx_save_model_state( mo, fname.c_str() );
   //mo.nc_save_model_state(fname.c_str());
  RETURNDIAG(rval=status);
}
// ............................................................... //
DIAGFUNC(diag_func_is_handle_valid) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=true);
  rval = is_value_valid_handle( getarg(0) );
  RETURNDIAG(rval);
}
// ............................................................... //
DIAGFUNC(diag_func_ncwritedims) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  // write the correct dimensions for the current model run
  //do nothing if model not compiled
  const int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  NcFile * ncfile = handler.get_nc_file(getarg(0));
  if(!ncfile) RETURNDIAG(rval=false);
  //NcDim * tindex = ncfile->add_dim("time");
  //NcDim * freeindex = 
  ncfile->add_dim("time");
  // NcDim * findex = 
  ncfile->add_dim("level", 8);
  //NcDim * fieldi =
  ncfile->add_dim("fieldi", numfields);
  //NcDim * fieldj = 
  ncfile->add_dim("fieldj", numfields);
  int ksize = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData().size(); 
  int nksize = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData().size(); 
  // NcDim * kindex = 
  ncfile->add_dim("kindex", ksize);
  // NcDim * nkindex= 
  ncfile->add_dim("nkindex", nksize);
  RETURNDIAG(rval=true);
};
// ............................................................... //
bool diag_func_ncwriteField(model_class & mo, psystem & p, 
    const diag_func_class & func_class, // needed to support __thisdiagfunc__ macro
    NcFile * ncfile, const qgcomplex * FIELD, const char * var1, const char * var2,
    int t1, int nctindex)
{
  const int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
  NcDim * fieldi = ncfile->get_dim("fieldi");
  NcDim * fieldj = ncfile->get_dim("fieldj");
  NcDim * kindex = ncfile->get_dim("kindex");
  NcDim * nkindex= ncfile->get_dim("nkindex");
  NcDim * time   = ncfile->get_dim("time");
  if(! (fieldi->is_valid() && fieldj->is_valid() && kindex->is_valid() && nkindex->is_valid() && time->is_valid() )){
    // didn't call adddims func
    set_mp_errmsg("Error in function '%s'. Necessary dimensions not found in file.",__thisdiagfunc__);
    return false;
  }
  const NcDim ** dimlist = new const NcDim*[5];
  dimlist[1] = fieldi; dimlist[2] = fieldj; dimlist[3] = kindex; dimlist[4] =nkindex;
  dimlist[0] = time;
  NcVar * rPhivar = ncgetaddvar(ncfile, var1, ncREAL, 5, dimlist);
  NcVar * iPhivar = ncgetaddvar(ncfile, var2, ncREAL, 5, dimlist);
  if(rPhivar == NULL || iPhivar == NULL){
    set_mp_errmsg("Error in function '%s'. Variable exists with inconsistent dimensions.",__thisdiagfunc__);
    return false;
  }
  if(nctindex < 0 || nctindex >= time->size() ){
    set_mp_errmsg("Error in function '%s'. Attempted to write data outside the temporal dimension size.",
        __thisdiagfunc__);
    return false;
  }
  //allocate the arrays
  REAL * rphi = new REAL[numfields*numfields * kindex->size() * nkindex->size()];
  REAL * iphi = new REAL[numfields*numfields * kindex->size() * nkindex->size()];
  int count = 0;
  // uncertain of the order here
  // current order is correct . . . . . (checked with NCL)
  for(int i = 0; i < numfields; i++){
    for(int j = 0; j < numfields; j++){
      for(int ki = 0; ki < kindex->size(); ki++){
        for(int nki = 0; nki < nkindex->size(); nki++){
          int index = ((nki*kindex->size()+ki)*fieldi->size() + i)*fieldj->size() + j; 
          //rank4 index(nki, ki, i, j);
          rphi[count] = FIELD[index].real();
          iphi[count] = FIELD[index].imag();
          count++; 
        } 
      }
    }
  }
  rPhivar->set_cur(nctindex, 0, 0, 0, 0);	
  iPhivar->set_cur(nctindex, 0, 0, 0, 0);
  rPhivar->put(rphi, 1, fieldi->size(), fieldj->size(), kindex->size(), nkindex->size());
  iPhivar->put(iphi, 1, fieldi->size(), fieldj->size(), kindex->size(), nkindex->size());
  delete [] rphi;
  delete [] iphi; 
  delete [] dimlist; 
  return true;
};
// ............................................................... //
DIAGFUNC(diag_func_ncgetmaxtime) {
  float_type rval;
  RETURNIFNOTCOMPILED(rval = 0.0);
  NcFile * ncfile = handler.get_nc_file(getarg(0));
  if(ncfile){
    const NcDim * time   = ncfile->get_dim("time");
    RETURNDIAG(rval=static_cast<float_type>( time->size() - 1 ));
  }else{
    mpError("Error in function '%s'. Invalid netcdf file handle supplied.",__thisdiagfunc__);
    RETURNDIAG(rval=-1);
  }
  RETURNDIAG(rval=-1); 
}
// ............................................................... //
DIAGFUNC(diag_func_ncaddtime) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  NcFile * ncfile = handler.get_nc_file(getarg(0));
  float_type tvalue = getarg(1).GetFloat(); 
  if(ncfile){
    const NcDim * time = ncfile->get_dim("time");
    NcVar * timevar = ncgetaddvar(ncfile, "time", ncREAL, 1, &(time));
    timevar->set_cur(time->size());
    REAL tval = static_cast<REAL>(tvalue);
    timevar ->put(&tval, 1, 0, 0, 0, 0);
    RETURNDIAG(static_cast<double>(time->size()));
  }else{
    RETURNDIAG(rval=false);
  }
  RETURNDIAG(rval=true); 
};
// ............................................................... //
DIAGFUNC(diag_func_ncwritePhi){ // this second one is for the netcdf located
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  NcFile * ncfile = handler.get_nc_file(getarg(0));
  if(ncfile){
    ncfile->sync();
    // check times
    psystem & p = mo.getPsystem();
    if( p.get_tindex() <= 0){
      set_mp_errmsg( "Model has not integrated forward." ); RETURNDIAG(rval=false);
    }
    int t1 = static_cast<int>(getarg(1).GetInteger()); // static_cast<int>(time);
    int t2 = static_cast<int>(getarg(2).GetInteger()); // static_cast<int>(time);
    if(t1 > p.get_tindex() || t1 < 0 || t2 > p.get_tindex() || t2 < 0){
      set_mp_errmsg( "Bad time argument in function '%s'.",__thisdiagfunc__);
      RETURNDIAG(rval=false);
    }
    const qgcomplex * cphi = (mo.getSsystem().get_phi_set(t1, t1).get_data());
    // write the current time to the time var
    RETURNDIAG(rval=diag_func_ncwriteField(mo, p, func_class, ncfile, cphi, "rPhi", "iPhi", t1, t2)); 
  }else{
    set_mp_errmsg( "Bad nc file handle passed to function %s.", __thisdiagfunc__);
    RETURNDIAG(rval=false);
  } 
};
// ............................................................... //
DIAGFUNC(diag_func_ncwriteG) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  NcFile * ncfile = handler.get_nc_file(getarg(0));
  if(ncfile){
    ncfile->sync();
    // check times
    psystem & p = mo.getPsystem();
    if( p.get_tindex() <= 0){
      set_mp_errmsg( "Model has not integrated forward." ); RETURNDIAG(rval=false);
    }
    int t1 = static_cast<int>(getarg(1).GetInteger()); // static_cast<int>(time);
    int t2 = static_cast<int>(getarg(2).GetInteger()); // static_cast<int>(time);
    if(t1 > p.get_tindex() || t1 < 0 || t2 > p.get_tindex() || t2 < 0){
      set_mp_errmsg( "Bad time argument in function '%s'.",__thisdiagfunc__);
      RETURNDIAG(rval=false);
    }
    const qgcomplex *  cG = (mo.getSsystem().get_phi_set(t1, t1).get_data());
    // write the current time to the time var
    RETURNDIAG(rval=diag_func_ncwriteField(mo, p, func_class, ncfile, cG, "rG", "iG", t1, t2)); 
  }else{
    set_mp_errmsg( "Bad nc file handle passed to function %s.", __thisdiagfunc__);
    RETURNDIAG(rval=false);
  } 
};
// ............................................................... //
DIAGFUNC(diag_func_ncclose) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  RETURNDIAG( rval =  handler.ncfile_close( getarg(0) ) );
}
// ............................................................... //
DIAGFUNC(diag_func_eval) { // evaluate contents of string
  RETURNIFNOTCOMPILED( get_default_value_of_type('v') );
  // evaluate the contents of a string
  const string & evalstr = getarg(0).GetString();
  ParserX p(pckNONE);
  mo.assignparserdata(p, "v", false);

  token_list expr;
  extract_words(evalstr, expr);

  // does not assign diag info
  p.SetExpr(expr);
  RETURNDIAG(p.Eval());
};
// ............................................................... //
DIAGFUNC(diag_func_readenv) {
  string rvalue;
  RETURNIFNOTCOMPILED(rvalue);
  const char * envstr = getenv(getarg(0).GetString().c_str());
  if(envstr){
    RETURNDIAG( rvalue=string(envstr) );
  }
  RETURNDIAG(rvalue);
};
// ............................................................... //
DIAGFUNC(diag_func_print_model_profiles) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  const qstream * fout = handler.get_file(getarg(0));
  if( ! fout ){
    RETURNDIAG(rval=false);
  }
  printMOProfiles(mo, *fout );
  RETURNDIAG(rval=true);
}
// ............................................................... //
DIAGFUNC(diag_func_getprofilemaxsize){
  RETURNDIAG( static_cast<int_type>( std::numeric_limits< INDEX >::max() ) );
}
// ............................................................... //
DIAGFUNC(diag_func_abort) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=false);
  exit(static_cast<int>(getarg(0).GetInteger()));
  RETURNDIAG(rval=true);
};
// ............................................................... //
DIAGFUNC(diag_func_end) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=true);
  raise(SIGTERM);
  RETURNDIAG(rval=true);
};
// ............................................................... //
DIAGFUNC(diag_func_geterrmsg){
  string_type rval;
  RETURNIFNOTCOMPILED(rval);
  RETURNDIAG( string_type( get_mp_errmsg() ) );
};
// ............................................................... //
DIAGFUNC(diag_func_print_model_macros) {
  bool_type rval;
  RETURNIFNOTCOMPILED(rval=true);
  const qstream * filehandle = handler.get_file(getarg(0));
  if( ! filehandle ) RETURNDIAG(rval=false);
  printMOMacros(mo, *filehandle );
  RETURNDIAG(rval=true);
}

#endif
